Optimizacion Heuristica

Integrantes

Juan Pablo Buitrago.

Jorge Andrés Higuita Monsalve.

Daniel Arango.

Parte 1: Optimizacion Numerica.

pip install deap

Requirement already satisfied: deap in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (1.3.3)Note: you may need to restart the kernel to use updated packages.

Requirement already satisfied: numpy in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from deap) (1.22.3)

pip install pyswarm

Requirement already satisfied: pyswarm in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (0.6)
Requirement already satisfied: numpy in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pyswarm) (1.22.3)
Note: you may need to restart the kernel to use updated packages.

pip install geopandas

Requirement already satisfied: geopandas in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (0.12.2)
Requirement already satisfied: shapely>=1.7 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (1.8.0)
Requirement already satisfied: pandas>=1.0.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (1.4.2)
Requirement already satisfied: fiona>=1.8 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (1.8.20)
Requirement already satisfied: pyproj>=2.6.1.post1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (3.1.0)
Requirement already satisfied: packaging in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from geopandas) (21.3)
Requirement already satisfied: attrs>=17 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (21.4.0)
Requirement already satisfied: certifi in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (2022.12.7)
Requirement already satisfied: click>=4.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (7.1.2)
Requirement already satisfied: cligj>=0.5 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (0.7.2)
Requirement already satisfied: click-plugins>=1.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (1.1.1)
Requirement already satisfied: six>=1.7 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (1.16.0)
Requirement already satisfied: munch in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (2.5.0)
Requirement already satisfied: setuptools in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from fiona>=1.8->geopandas) (61.2.0)
Requirement already satisfied: python-dateutil>=2.8.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pandas>=1.0.0->geopandas) (2.8.2)
Requirement already satisfied: numpy>=1.18.5 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pandas>=1.0.0->geopandas) (1.22.3)
Requirement already satisfied: pytz>=2020.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pandas>=1.0.0->geopandas) (2022.1)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from packaging->geopandas) (3.0.9)
Note: you may need to restart the kernel to use updated packages.

pip install pygad

Requirement already satisfied: pygad in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (2.19.2)
Requirement already satisfied: numpy in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pygad) (1.22.3)
Requirement already satisfied: matplotlib in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pygad) (3.5.2)
Requirement already satisfied: cloudpickle in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from pygad) (2.0.0)
Requirement already satisfied: fonttools>=4.22.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (4.33.3)
Requirement already satisfied: pillow>=6.2.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (9.0.1)
Requirement already satisfied: packaging>=20.0 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (21.3)
Requirement already satisfied: python-dateutil>=2.7 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (2.8.2)
Requirement already satisfied: kiwisolver>=1.0.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (1.4.3)
Requirement already satisfied: pyparsing>=2.2.1 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (3.0.9)
Requirement already satisfied: cycler>=0.10 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from matplotlib->pygad) (0.11.0)
Requirement already satisfied: six>=1.5 in c:\users\hp\anaconda3\envs\tf_gpu\lib\site-packages (from python-dateutil>=2.7->matplotlib->pygad) (1.16.0)
Note: you may need to restart the kernel to use updated packages.

#Importamos las librerias Necesarias
import random
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from matplotlib.colors import LogNorm
import matplotlib.animation as animation
from IPython.display import HTML
from deap import algorithms, base, creator, tools
from pyswarm import pso
from scipy.optimize import minimize
import geopandas as gpd
from shapely.geometry import Point
import pyproj
import pygad 
from pygad import gann
from matplotlib import animation, rc
import pandas as pd
import itertools
import imageio
import os
from matplotlib.animation import FuncAnimation

Función Rastrigin

# define la función de Rastrigin de dos dimensiones
def f_rastrigin2d(x1, x2):
  y = 20 + x1 ** 2 - 10 * np.cos(2 * np.pi * x1) + x2 ** 2 - 10 * np.cos(2 * np.pi * x2)
  return y

# define la función de Rastrigin de tres dimensiones
def f_rastrigin3d(x1,x2,x3):
  y = 20 + x1 ** 2 - 10 * np.cos(2 * np.pi * x1) + x2 ** 2 - 10 * np.cos(2 * np.pi * x2) + x3 ** 2 - 10 * np.cos(2 * np.pi * x3)
  return y

n = 50  # define el número de puntos en la malla

# crea una lista de valores igualmente espaciados desde -5.12 a 5.12
x1 = np.linspace(-5.12, 5.12, n)
x2 = np.linspace(-5.12, 5.12, n)

# crea una malla de valores de entrada 2D utilizando los valores de x1 y x2
X1, X2 = np.meshgrid(x1, x2)

# aplana la malla y concatena en una matriz de entrada X
X = np.vstack((X1.flatten(), X2.flatten())).T

# calcula el valor de la función de Rastrigin de dos dimensiones para cada punto en la matriz X
z = f_rastrigin2d(X[:, 0], X[:, 1])

# reformatea el resultado en una matriz 2D de tamaño n x n
Z = z.reshape(n, n)

Funcion Rastrigin 2D

# Crear un gráfico de contorno con valores numéricos
levels = np.arange(0, 90, 10)
cp = plt.contourf(x1, x2, Z, levels=levels, cmap='plasma')
plt.colorbar(cp)

# Agregar líneas de contorno
plt.contour(x1, x2, Z, levels=levels, colors='k', linewidths=0.5)

# Título y etiquetas de los ejes
plt.title("Función de Rastrigin")
plt.xlabel("x1")
plt.ylabel("x2")

# Mostrar el gráfico
plt.show()

Función Rastrigin 3D

# Cálculo de la función
z = f_rastrigin2d(X1, X2)

# Creación de la figura
fig = go.Figure(data=[go.Surface(x=X1, y=X2, z=z)])

# Personalización de la figura
fig.update_layout(
    title='Función de Rastrigin 2D',
    scene = dict(
        xaxis_title='X1',
        yaxis_title='X2',
        zaxis_title='Y',
        aspectmode='cube',
        camera=dict(
            eye=dict(x=1.75, y=-1.75, z=1.25)
        )
    )
)

# Mostrar la figura
fig.show()

Optimización funcion Rastrigin

Descenso por gradiente

En 2 dimensiones

Vamos a usar el metodo de descenso de gradiente llamado metodo de Newton- Raphson

# Define una función que calcula la función de Rastrigin dada una entrada x
def rastrigin(x):
    n = len(x) # Obtiene la longitud de la entrada x
    return 10*n + np.sum(x**2 - 10*np.cos(2*np.pi*x)) # Calcula y devuelve el valor de la función de Rastrigin

# Define una función que calcula el gradiente de la función de Rastrigin dada una entrada x
def gradient(x):
    n = len(x) # Obtiene la longitud de la entrada x
    grad = np.zeros(n) # Crea un vector gradiente de ceros de longitud n
    for i in range(n): # Itera sobre cada entrada de la entrada x
        grad[i] = 2*x[i] + 20*np.pi*np.sin(2*np.pi*x[i]) # Calcula la entrada correspondiente en el vector gradiente
    return grad # Devuelve el vector gradiente

# Define una función que calcula la matriz Hessiana de la función de Rastrigin dada una entrada x
def hessian(x):
    n = len(x) # Obtiene la longitud de la entrada x
    hess = np.zeros((n, n)) # Crea una matriz Hessiana de ceros de tamaño n x n
    for i in range(n): # Itera sobre cada entrada de la entrada x
        for j in range(n): # Itera sobre cada entrada de la entrada x
            if i == j: # Si i y j son iguales
                hess[i][i] = 2 + 40*np.pi**2*np.cos(2*np.pi*x[i]) # Calcula la entrada correspondiente en la matriz Hessiana
            else:
                hess[i][j] = 0 # De lo contrario, establece la entrada correspondiente en la matriz Hessiana a cero
    return hess # Devuelve la matriz Hessiana

# Define una función que optimiza la función de Rastrigin dada una conjetura inicial, un número máximo de iteraciones y una tolerancia
def optimize_rastrigin(initial_guess, max_iterations, tolerance=1e-6):
    x = initial_guess # Establece la conjetura inicial como el valor actual de x
    path = [x] # Crea una lista que almacena los valores de x en cada iteración
    for i in range(max_iterations): # Itera sobre el número máximo de iteraciones
        grad = gradient(x) # Calcula el gradiente de la función de Rastrigin en el valor actual de x
        hess = hessian(x) # Calcula la matriz Hessiana de la función de Rastrigin en el valor actual de x
        step = np.linalg.solve(hess, -grad) # Calcula el paso de optimización utilizando la solución de la ecuación Hessiana
        x_new = x + step # Calcula el nuevo valor de x utilizando el paso de optimización
        if np.linalg.norm(x_new - x) < tolerance: # Si la norma de la diferencia entre el nuevo y el valor actual de x es menor que la tolerancia, la optimización se ha completado
            break
        x = x_new # De lo contrario, establece el nuevo valor de x como el valor actual de x y continúa iterando
        path.append(x) # Agrega el nuevo valor de x a la lista de valores de x en cada iteración
    return x, path # Devuelve el valor final

rastrigin_init_point_2D = np.random.uniform(-5.12, 5.12, size=2) #Se genera un número aleatorio y se acota entre el dominio de la función
max_iterations = 1000  
rastrigin_min_2D,rastrigin_path_2D = optimize_rastrigin(rastrigin_init_point_2D, 1000)

#Impresión de resultados
print(f"Init value: {rastrigin_init_point_2D}")
print(f"Minimum: {rastrigin_min_2D}")
print(f"Function value at minimum: {rastrigin(rastrigin_min_2D)}")
print("Optimization path: ")
for point in rastrigin_path_2D:
    print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {rastrigin(point)}")

    # Generate a grid of x, y values
rtg_x_min, rtg_x_max, rtg_x_maxvalues = -5.12, 5.12, 100
rtg_y_min, rtg_y_max, rtg_y_maxvalues = -5.12, 5.12, 100

x_rastrigin = np.linspace(rtg_x_min, rtg_x_max, rtg_x_maxvalues)
y_rastrigin = np.linspace(rtg_y_min, rtg_y_max, rtg_y_maxvalues)
X_rastrigin, Y_rastrigin = np.meshgrid(x_rastrigin, y_rastrigin)
# Compute the Rastrigin function for each (x, y) pair
Z_rastrigin = np.zeros_like(X_rastrigin)
for i in range(X_rastrigin.shape[0]):
    for j in range(X_rastrigin.shape[1]):
        Z_rastrigin[i, j] = rastrigin(np.array([X_rastrigin[i, j], Y_rastrigin[i, j]]))

rastrigin_minimum= rastrigin_min_2D.reshape(-1, 1)

#Test plot

rastrigin_surface_fig = plt.figure(figsize=(10,8))
rastrigin_surface_axes = plt.axes(projection='3d', elev=50, azim=-50)

rastrigin_surface_axes.plot_surface(X_rastrigin, Y_rastrigin, Z_rastrigin, norm=LogNorm(), rstride=1, cstride=1, 
                edgecolor='none', alpha=.2, cmap=plt.cm.jet)
rastrigin_surface_axes.plot(*rastrigin_minimum, rastrigin(rastrigin_minimum), 'b*', markersize=10)

rastrigin_surface_axes.set_xlabel('$x$')
rastrigin_surface_axes.set_ylabel('$y$')
rastrigin_surface_axes.set_zlabel('$z$')

rastrigin_surface_axes.set_xlim((rtg_x_min, rtg_x_max))
rastrigin_surface_axes.set_ylim((rtg_y_min, rtg_y_max))

plt.show()

# Create a 2D contour plot of the Rastrigin function and plot the ptimization path
rastrigin_contour_fig, rastrigin_contour_axes = plt.subplots(figsize=(10,8))
contour = rastrigin_contour_axes.contour(X_rastrigin, Y_rastrigin, Z_rastrigin, levels=20, cmap='coolwarm')
rastrigin_contour_axes.set_xlabel('x')
rastrigin_contour_axes.set_ylabel('y')

for point in rastrigin_path_2D:
  rastrigin_contour_axes.plot(*point, 'o-',color='blue',alpha=1)  

# Add a marker for the global minimum at (0, 0)
rastrigin_contour_axes.scatter(*rastrigin_minimum, marker='*', s=500, color='red')
plt.colorbar(contour)
plt.show()

Init value: [-0.28367766 -1.40353059]
Minimum: [-0.99495911 -1.50764073]
Function value at minimum: 23.256417941943447
Optimization path: 
Punto evaluado:[-0.28367766 -1.40353059]-- Valor de la función en el punto: 32.369204323421876
Punto evaluado:[-1.04988166 -1.5232434 ]-- Valor de la función en el punto: 23.803209023269208
Punto evaluado:[-0.99301162 -1.50755286]-- Valor de la función en el punto: 23.25716802783503
Punto evaluado:[-0.99495911 -1.50764073]-- Valor de la función en el punto: 23.256417941943447

Animacion descenso

# Define the function to update the animation
def update(frame):
    x = rastrigin_path_2D[frame]
    line.set_data(x[0], x[1])
    if frame > 0:
        line_trace.set_data([rastrigin_path_2D[i][0] for i in range(frame+1)], [rastrigin_path_2D[i][1] for i in range(frame+1)])
    return line, line_trace
# Plot the contour and the starting point
rastrigin_contour_fig, rastrigin_contour_axes = plt.subplots(figsize=(10,8))
rastrigin_contour_axes.set_title('Optimization Path Animation')
rastrigin_contour_axes.set_xlabel('x')
rastrigin_contour_axes.set_ylabel('y')
rastrigin_contour = rastrigin_contour_axes.contour(X_rastrigin, Y_rastrigin, Z_rastrigin, levels=20, cmap='coolwarm')
rastrigin_contour_axes.plot(rastrigin_init_point_2D[0], rastrigin_init_point_2D[1], 'o', color='blue')

# Set up the animation

line, = rastrigin_contour_axes.plot([], [], 'o-', color='red', lw=10)
line_trace, = rastrigin_contour_axes.plot([], [], '-', color='green', lw=4)
ani = animation.FuncAnimation(rastrigin_contour_fig, update, frames=len(rastrigin_path_2D), interval=500, blit=True)

HTML(ani.to_jshtml())

Once Loop Reflect

En 3 dimensiones

rastrigin_init_point_3D = np.random.uniform(-5.12, 5.12, size=3) #Se genera un número aleatorio y se acota entre el dominio de la función
rastrigin_min_3D,rastrigin_path_3D = optimize_rastrigin(rastrigin_init_point_3D, 1000)

#Impresión de resultados
print("Init value", rastrigin_init_point_3D)
print("Minimum:", rastrigin_min_3D)
print("Function value at minimum:", rastrigin(rastrigin_min_3D))
print("Optimization path: ")
for point in rastrigin_path_3D:
    print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {rastrigin(point)}")

Init value [-2.87663391  2.66267468  0.51743905]
Minimum: [-2.9848557   2.51274332  0.50254604]
Function value at minimum: 55.487715409989356
Optimization path: 
Punto evaluado:[-2.87663391  2.66267468  0.51743905]-- Valor de la función en el punto: 53.64496918876429
Punto evaluado:[-3.01119854  2.42587958  0.50249088]-- Valor de la función en el punto: 54.16324400210692
Punto evaluado:[-2.98482268  2.52015726  0.50254604]-- Valor de la función en el punto: 55.47697082306979
Punto evaluado:[-2.9848557   2.51272392  0.50254604]-- Valor de la función en el punto: 55.487715336295715
Punto evaluado:[-2.9848557   2.51274332  0.50254604]-- Valor de la función en el punto: 55.487715409989356

Algoritmos evolutivos

En 2 dimensiones

# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0

# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5

# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)

# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initCycle, creator.Individuo,
                 (toolbox.attr_float, toolbox.attr_float), n=1)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)

# Definir la función de evaluación
def evaluar(individuo):
    return f_rastrigin2d(individuo[0], individuo[1]),

# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)

# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)

# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
    # Seleccionar a los padres
    padres = toolbox.select(poblacion, len(poblacion))
    # Clonar los padres para aplicar operadores genéticos
    hijos = list(map(toolbox.clone, padres))

    # Aplicar cruce y mutación a los hijos
    for hijo in hijos:
        if random.random() < probabilidad_cruce:
            toolbox.mate(hijo, random.choice(hijos))
            del hijo.fitness.values

        if random.random() < probabilidad_mutacion:
            toolbox.mutate(hijo)
            del hijo.fitness.values

    # Evaluar la aptitud de los nuevos individuos
    individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
    aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
    for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
        individuo.fitness.values = aptitud

    # Reemplazar la población anterior por la nueva generación
    poblacion = hijos

# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])

Mejor solución encontrada: [8.604753801223136e-09, 3.210492889238655e-09]
Valor de la función objetivo: 1.5987211554602254e-14

En 3 dimensiones

# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0

# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5

# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)

# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initRepeat, creator.Individuo,
                 toolbox.attr_float, n=3)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)

# Definir la función de evaluación
def evaluar(individuo):
    return f_rastrigin3d(individuo[0], individuo[1], individuo[2]),

# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)

# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)

# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
    # Seleccionar a los padres
    padres = toolbox.select(poblacion, len(poblacion))
    # Clonar los padres para aplicar operadores genéticos
    hijos = list(map(toolbox.clone, padres))

    # Aplicar cruce y mutación a los hijos
    for hijo in hijos:
        if random.random() < probabilidad_cruce:
            toolbox.mate(hijo, random.choice(hijos))
            del hijo.fitness.values

        if random.random() < probabilidad_mutacion:
            toolbox.mutate(hijo)
            del hijo.fitness.values

    # Evaluar la aptitud de los nuevos individuos
    individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
    aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
    for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
        individuo.fitness.values = aptitud

    # Reemplazar la población anterior por la nueva generación
    poblacion = hijos

# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])

Mejor solución encontrada: [0.019417809579997895, 0.009269224741063526, 0.0009761931431016592]
Valor de la función objetivo: -9.908058456113144

c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning:

A class named 'FitnessMin' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.

c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning:

A class named 'Individuo' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.

Optimizacion de particulas

En 2 dimensiones

#Definimos la funcion rastringin para los calculos
def f(x, y):
    return (20 + x ** 2 - 10 * np.cos(2 * np.pi * x) + y ** 2 - 10 * np.cos(2 * np.pi * y))


#Calculamos los datos X Y de la funcion
x, y = np.array(np.meshgrid(np.linspace(-5.12, 5.12, 50),
                np.linspace(-5.12, 5.12, 50)))
#Evalamos los datos X Y
z = f(x, y)

#Encontramos el minimo Global de la funcion
x_min = x.ravel()[z.argmin()]
y_min = y.ravel()[z.argmin()]
# Se inicializan los Hiperparametros del algoritmo PSO
c1 = c2 = 0.1  # Factores de Acelearacion
w = 0.8  # Coeficiente de Inercia

# Creamos las particulas
n_particles = 35
np.random.seed(100)
X = np.random.rand(2, n_particles) * 5  # Posiscion Inicial de las particulas
V = np.random.randn(2, n_particles) * 0.5  # velocidad de las particulas

# Inicializacion de datos
pbest = X  # Mejor Posicion de la particula
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()]  # Mejor Posicion del Enjambre
gbest_obj = pbest_obj.min()
#funciond de actualizacion de los puntos del enjambre


def update():
    global V, X, pbest, pbest_obj, gbest, gbest_obj
    # actualizacion de parametros
    r1, r2 = np.random.rand(2)
    # Actualizacion de velocidades
    V = w * V + c1*r1*(pbest - X) + c2*r2*(gbest.reshape(-1, 1)-X)
    X = X + V  # Actualizacion de Posicion de particulas
    obj = f(X[0], X[1])  # Nueva Posicion
    # Se evalua si la posicion nueva 'obj' es mejor a la anterior
    pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]
    pbest_obj = np.array([pbest_obj, obj]).min(axis=0)
    gbest = pbest[:, pbest_obj.argmin()]
    gbest_obj = pbest_obj.min()


# Visualizacion del grafico de Grienwank Y los puntos del enjambre
fig, ax = plt.subplots(figsize=(16, 8))
fig.set_tight_layout(True)
img = ax.imshow(z, extent=[-5.12, 5.12, -5.12, 5.12],
                origin='lower', cmap='plasma', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([x_min], [y_min], marker='x', markersize=10,
        color="white")  # Minimo Global
contours = ax.contour(x, y, z, 10, colors='black', alpha=0.4)
ax.clabel(contours, inline=True, fontsize=8, fmt="%.0f")
# mejor posision del punto
pbest_plot = ax.scatter(
    pbest[0], pbest[1], marker='o', color='black', alpha=0.5)
p_plot = ax.scatter(X[0], X[1], marker='o', color='green', alpha=0.5)  # Puntos
p_arrow = ax.quiver(X[0], X[1], V[0], V[1], color='green', width=0.005,
                    angles='xy', scale_units='xy', scale=1)  # flecha de los puntos
gbest_plot = plt.scatter([gbest[0]], [gbest[1]], marker='*',
                         s=100, color='red', alpha=0.4)  # Mejor posicion del enjambre
ax.set_xlim([-5.12, 5.12])
ax.set_ylim([-5.12, 5.12])


def animate(i):
    # Animacion de los puntos en moviento
    title = 'Iteration {:02d}'.format(i)
    # Update params
    update()
    # Set picture
    ax.set_title(title)
    pbest_plot.set_offsets(pbest.T)
    p_plot.set_offsets(X.T)
    p_arrow.set_offsets(X.T)
    p_arrow.set_UVC(V[0], V[1])
    gbest_plot.set_offsets(gbest.reshape(1, -1))
    return ax, pbest_plot, p_plot, p_arrow, gbest_plot


anim = FuncAnimation(fig, animate, frames=list(
    range(1, 50)), interval=500, blit=False, repeat=True)
# Se guarda la animacion en formato gif
anim.save("PSORAS.gif", dpi=120, writer="imagemagick")
# Informacion genreal del algoritmo luego de optimizar
print("PSO found best solution at f({})={}".format(gbest, gbest_obj))
print("Global optimal at f({})={}".format([x_min,y_min], f(x_min,y_min)))

MovieWriter imagemagick unavailable; using Pillow instead.

PSO found best solution at f([-0.0005921  -0.00046606])=0.00011264565412716365
Global optimal at f([-0.9404081632653059, -0.9404081632653059])=3.1543848912857424

Gif Funcion Rastrigin

from IPython.display import Image, display

Image(filename='PSORAS.gif')

<IPython.core.display.Image object>

En 3 dimensiones

lb = [-10, -10, -10]  # límite inferior de las variables
ub = [10, 10, 10]     # límite superior de las variables

x_opt, f_opt = pso(lambda x: f_rastrigin3d(x[0], x[1], x[2]), lb, ub)

print("El punto óptimo es:", x_opt)

Stopping search: Swarm best objective change less than 1e-08
El punto óptimo es: [ 2.05480392e-05 -9.57955878e-06  1.71508643e-05]

Optimizacion por evolucion diferencial

En 2 dimensiones

# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones

# Definir el límite inferior y superior para los valores de x e y
lim_inf = -10
lim_sup = 10

# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
    x = random.uniform(lim_inf, lim_sup)
    y = random.uniform(lim_inf, lim_sup)
    pop.append((x, y))

# Realizar la evolución diferencial
for gen in range(ngen):
    for i in range(NP):
        # Seleccionar 3 individuos de forma aleatoria
        a, b, c = random.sample(pop, 3)

        # Seleccionar un número aleatorio para cada variable
        R1, R2 = random.random(), random.random()

        # Realizar la mutación diferencial
        x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
        y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)

        # Realizar la recombinación
        if random.random() < CR:
            x = x_new
            y = y_new
        else:
            x = a[0]
            y = a[1]

        # Evaluar la función en el nuevo individuo
        fitness_new = f_rastrigin2d(x, y)

        # Reemplazar el individuo anterior si el nuevo individuo es mejor
        if fitness_new < f_rastrigin2d(pop[i][0], pop[i][1]):
            pop[i] = (x, y)

# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_rastrigin2d(x[0], x[1]))
print('Valor mínimo de la función: ', f_rastrigin2d(best_ind[0], best_ind[1]))
print('Variables que producen el valor mínimo: ', best_ind)

Valor mínimo de la función:  0.0
Variables que producen el valor mínimo:  (-1.6130253689405727e-09, -1.3999432868862866e-09)

En 3 dimensiones

# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones

# Definir el límite inferior y superior para los valores de x, y, y z
lim_inf = -10
lim_sup = 10

# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
    x = random.uniform(lim_inf, lim_sup)
    y = random.uniform(lim_inf, lim_sup)
    z = random.uniform(lim_inf, lim_sup)
    pop.append((x, y, z))

# Realizar la evolución diferencial
for gen in range(ngen):
    for i in range(NP):
        # Seleccionar 3 individuos de forma aleatoria
        a, b, c = random.sample(pop, 3)

        # Seleccionar un número aleatorio para cada variable
        R1, R2, R3 = random.random(), random.random(), random.random()

        # Realizar la mutación diferencial
        x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
        y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)
        z_new = np.clip(a[2] + F*(b[2] - c[2]), lim_inf, lim_sup)

        # Realizar la recombinación
        if random.random() < CR:
            x = x_new
            y = y_new
            z = z_new
        else:
            x = a[0]
            y = a[1]
            z = a[2]

        # Evaluar la función en el nuevo individuo
        fitness_new = f_rastrigin3d(x, y, z)

        # Reemplazar el individuo anterior si el nuevo individuo es mejor
        if fitness_new < f_rastrigin3d(pop[i][0], pop[i][1], pop[i][2]):
            pop[i] = (x, y, z)

# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_rastrigin3d(x[0], x[1], x[2]))
print('Valor mínimo de la función: ', f_rastrigin3d(best_ind[0], best_ind[1], best_ind[2]))
print('Variables que producen el valor mínimo: ', best_ind)

Valor mínimo de la función:  -5.2820006782696485
Variables que producen el valor mínimo:  (0.04069343900956812, 0.09806400831968248, 1.0842032194549434)
 -5.2820006782696485
Variables que producen el valor mínimo:  (0.04069343900956812, 0.09806400831968248, 1.0842032194549434)

Funcion Griewank

def f_griewank2d(x, y):
    return 1 + (x**2 + y**2)/4000 - np.cos(x)*np.cos(y/np.sqrt(2))

def f_griewank3d(x, y, z):
    return 1 + (x**2 + y**2 + z**2)/4000 - np.cos(x)*np.cos(y/np.sqrt(2))*np.cos(z/np.sqrt(3))

n = 50

x1 = np.linspace(-600, 600, n)
x2 = np.linspace(-600, 600, n)

X1, X2 = np.meshgrid(x1, x2)
X = np.vstack((X1.flatten(), X2.flatten())).T
z = f_griewank2d(X[:, 0], X[:, 1])
Z = z.reshape(n, n)

Funcion Griewank 2D

# Crear un gráfico de contorno con valores numéricos
levels = np.arange(0, 201, 10)
cp = plt.contourf(x1, x2, Z, levels=levels, cmap='plasma')
plt.colorbar(cp)

# Agregar líneas de contorno
plt.contour(x1, x2, Z, levels=levels, colors='k', linewidths=0.5)

# Título y etiquetas de los ejes
plt.title("Función de Griewank")
plt.xlabel("x1")
plt.ylabel("x2")

# Mostrar el gráfico
plt.show()

Funcion Griewank 3D

# Cálculo de la función
z = f_griewank2d(X1, X2)

# Creación de la figura
fig = go.Figure(data=[go.Surface(x=X1, y=X2, z=z)])

# Personalización de la figura
fig.update_layout(
    title='Función de Griewank 2D',
    scene = dict(
        xaxis_title='X1',
        yaxis_title='X2',
        zaxis_title='griewank(X1,X2)',
        aspectmode='cube',
        camera=dict(
            eye=dict(x=1.75, y=-1.75, z=1.25)
        )
    )
)

# Mostrar la figura
fig.show()

Optimización

Descenso por gradiente

En 2 dimensiones

def griewank(x):
    # Esta función calcula el valor de la función Griewank en un punto dado x
    n = x.shape[0]
    sum = np.sum(x**2)
    prod = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1))))
    return 1 + sum / 4000 - prod

def griewank_gradient(x):
    # Esta función calcula el gradiente de la función Griewank en un punto dado x
    n = x.shape[0]
    prod = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1))))
    grad = np.zeros(n)
    for i in range(n):
        grad[i] = x[i] / 2000 + (1 / np.sqrt(i+1)) * prod * np.sin(x[i] / np.sqrt(i+1))
    return grad

def griewank_hessian(x):
    # Esta función calcula la matriz hessiana de la función Griewank en un punto dado x
    n = x.shape[0]
    prod = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1))))
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(i, n):
            if i == j:
                H[i][j] = 1 / 2000 + (1 / (i+1)) * prod * (1 - (np.sin(x[i] / np.sqrt(i+1)))**2)
            else:
                pij = np.prod(np.cos(x / np.sqrt(np.arange(1, n+1)))[:i]) * np.prod(np.cos(x / np.sqrt(np.arange(1, n+1)))[i+1:j])
                H[i][j] = (-1 / ((i-j)**2)) * (np.sin(x[i] / np.sqrt(i+1))) * (np.sin(x[j] / np.sqrt(j+1))) * prod * pij
                H[j][i] = H[i][j]
    return H

def newton_raphson_optimization(x, max_iter=100, eps=1e-6):
    # Esta función aplica el método de Newton-Raphson para encontrar el mínimo de la función Griewank
    # a partir de un punto inicial x
    path = [x]
    for i in range(max_iter):
        grad = griewank_gradient(x)
        H = griewank_hessian(x)
        delta_x = np.linalg.solve(H, -grad)
        x = x + delta_x
        path.append(x)
        if np.linalg.norm(delta_x) < eps:
            break
    return x, path

griewank_scale = 600 / 3  # Scale factor for desired range of [-600, 600] (3 standard deviations)
griewank_init_point_2D = griewank_scale * np.random.randn(2)
griewank_init_point_2D = np.clip(griewank_init_point_2D, -600, 600)  # Clip values outside of the desired range
griewank_min_2D, griewank_path_2D = newton_raphson_optimization(griewank_init_point_2D, max_iter= 10000, )

#Impresión de resultados
print(f"Init value: {griewank_init_point_2D}")
print(f"Minimum: {griewank_min_2D}")
print(f"Function value at minimum: {griewank(griewank_min_2D)}")
print("Optimization path: ")
for point in griewank_path_2D:
    print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {griewank(point)}")

Init value: [-268.45690276 -328.4636416 ]
Minimum: [-219.79895285 -239.66954303]
Function value at minimum: 26.459604466919707
Optimization path: 
Punto evaluado:[-268.45690276 -328.4636416 ]-- Valor de la función en el punto: 46.1345318918537
Punto evaluado:[-236.8090539  -316.80549964]-- Valor de la función en el punto: 39.89838123232788
Punto evaluado:[-236.92365268 -316.33887784]-- Valor de la función en el punto: 39.838259647444545
Punto evaluado:[-237.09813642 -315.64651031]-- Valor de la función en el punto: 39.871017739133634
Punto evaluado:[-266.2349533  -320.68007028]-- Valor de la función en el punto: 45.019080016296016
Punto evaluado:[-258.46719463 -312.11084895]-- Valor de la función en el punto: 41.590905278617775
Punto evaluado:[-262.4683787  -317.27379128]-- Valor de la función en el punto: 43.427757027481704
Punto evaluado:[-257.43322075 -312.71416805]-- Valor de la función en el punto: 41.66878110225802
Punto evaluado:[-262.42626699 -282.19693939]-- Valor de la función en el punto: 38.120281633940586
Punto evaluado:[-290.61622545 -308.67320533]-- Valor de la función en el punto: 45.932803518253735
Punto evaluado:[-209.4212731  -234.30413624]-- Valor de la función en el punto: 25.36090527066114
Punto evaluado:[-211.75824503 -237.01649732]-- Valor de la función en el punto: 26.11865563006501
Punto evaluado:[-212.02304479 -236.8364216 ]-- Valor de la función en el punto: 26.241675864745016
Punto evaluado:[-219.94341084 -242.45003907]-- Valor de la función en el punto: 28.008635536129148
Punto evaluado:[-220.23055931 -236.56649934]-- Valor de la función en el punto: 27.795748849615975
Punto evaluado:[-219.213345   -240.05651962]-- Valor de la función en el punto: 26.65816286357996
Punto evaluado:[-220.12552779 -239.47943438]-- Valor de la función en el punto: 26.520355625595574
Punto evaluado:[-219.76322165 -239.72024788]-- Valor de la función en el punto: 26.460817235236306
Punto evaluado:[-219.79852135 -239.6708955 ]-- Valor de la función en el punto: 26.459602829915493
Punto evaluado:[-219.79894308 -239.66956946]-- Valor de la función en el punto: 26.45960442541122
Punto evaluado:[-219.79895265 -239.66954355]-- Valor de la función en el punto: 26.459604466104288
Punto evaluado:[-219.79895285 -239.66954303]-- Valor de la función en el punto: 26.459604466919707

# Generate a grid of x, y values
gwk_x_min, gwk_x_max, gwk_x_maxvalues = -600,600, 600
gwk_y_min, gwk_y_max, gwk_y_maxvalues = -600, 600, 600

x_griewank = np.linspace(gwk_x_min, gwk_x_max, gwk_x_maxvalues)
y_griewank = np.linspace(gwk_y_min, gwk_y_max, gwk_y_maxvalues)
X_griewank, Y_griewank = np.meshgrid(x_griewank, y_griewank)
# Compute the Rastrigin function for each (x, y) pair
Z_griewank = np.zeros_like(X_griewank)
for i in range(X_griewank.shape[0]):
    for j in range(X_griewank.shape[1]):
        Z_griewank[i, j] = griewank(np.array([X_griewank[i, j], Y_griewank[i, j]]))


griewank_minimum= griewank_min_2D.reshape(-1, 1)

#Plot Griewank surface and minimum point founded.

griewank_surface_fig = plt.figure(figsize=(10,8))
griewank_surface_axes = plt.axes(projection='3d', elev=50, azim=-50)

griewank_surface_axes.plot_surface(X_griewank, Y_griewank, Z_griewank, norm=LogNorm(), rstride=1, cstride=1, 
                edgecolor='none', alpha=1, cmap=plt.cm.jet)
griewank_surface_axes.plot(*griewank_minimum, griewank(griewank_minimum), 'b*', markersize=10)

griewank_surface_axes.set_xlabel('$x$')
griewank_surface_axes.set_ylabel('$y$')
griewank_surface_axes.set_zlabel('$z$')

griewank_surface_axes.set_xlim((gwk_x_min, gwk_x_max))
griewank_surface_axes.set_ylim((gwk_y_min, gwk_y_max))

plt.show()

# Create a 2D contour plot of the Griewank function and plot the optimization path
griewank_contour_fig, griewank_contour_axes = plt.subplots(figsize=(10,8))
griewank_contour = griewank_contour_axes.contour(X_griewank, Y_griewank, Z_griewank, levels=20, cmap='coolwarm')
griewank_contour_axes.set_xlabel('x')
griewank_contour_axes.set_ylabel('y')

for point in griewank_path_2D:
  griewank_contour_axes.plot(*point, 'o-',color='blue',alpha=1)  

# Add a marker for the global minimum at (0, 0)
griewank_contour_axes.scatter(*griewank_minimum, marker='*', s=500, color='red')
plt.colorbar(griewank_contour)
plt.show()

Animación descenso

# Define the function to update the animation
def update(frame):
    x = griewank_path_2D[frame]
    line.set_data(x[0], x[1])
    if frame > 0:
        line_trace.set_data([griewank_path_2D[i][0] for i in range(frame+1)], [griewank_path_2D[i][1] for i in range(frame+1)])
    return line, line_trace
# Set the starting point for the optimization
# x_start = np.array([-4, 4])

# Run the optimization
# x_opt, path = newton_raphson_optimization(x_start)

# Plot the contour and the starting point
griewank_contour_fig, griewank_contour_axes = plt.subplots(figsize=(10,8))
griewank_contour_axes.set_title('Optimization Path Animation')
griewank_contour_axes.set_xlabel('x')
griewank_contour_axes.set_ylabel('y')
griewank_contour = griewank_contour_axes.contour(X_griewank, Y_griewank, Z_griewank, levels=20, cmap='coolwarm')
griewank_contour_axes.plot(griewank_init_point_2D[0], griewank_init_point_2D[1], 'o', color='blue')

# Set up the animation

line, = griewank_contour_axes.plot([], [], 'o-', color='red', lw=2)
line_trace, = griewank_contour_axes.plot([], [], '-', color='green', lw=4)
ani = animation.FuncAnimation(griewank_contour_fig, update, frames=len(griewank_path_2D), interval=200, blit=True)

HTML(ani.to_jshtml())

Once Loop Reflect

En 3 dimensiones

griewank_scale = 600 / 3  # Scale factor for desired range of [-600, 600] (3 standard deviations)
griewank_init_point_3D = griewank_scale * np.random.randn(3)
griewank_init_point_3D = np.clip(griewank_init_point_3D, -600, 600)  # Clip values outside of the desired range
griewank_min_3D, griewank_path_3D = newton_raphson_optimization(griewank_init_point_3D, max_iter= 10000, )

#Impresión de resultados
print(f"Init value: {griewank_init_point_3D}")
print(f"Minimum: {griewank_min_3D}")
print(f"Function value at minimum: {griewank(griewank_min_3D)}")
print("Optimization path: ")
for point in griewank_path_3D:
    print(f"Punto evaluado:{point}--", f"Valor de la función en el punto: {griewank(point)}")

Init value: [ 21.2161394  465.38624143  20.40039014]
Minimum: [ 25.11933844 461.55675686  21.7308067 ]
Function value at minimum: 53.59734061619003
Optimization path: 
Punto evaluado:[ 21.2161394  465.38624143  20.40039014]-- Valor de la función en el punto: 55.007908731012144
Punto evaluado:[ 21.15121841 466.11706225  24.21638111]-- Valor de la función en el punto: 55.47495863734799
Punto evaluado:[ 20.70286608 463.18230165  17.23096299]-- Valor de la función en el punto: 54.646545826382294
Punto evaluado:[ 24.06364336 465.14289142  17.82271096]-- Valor de la función en el punto: 55.135178039985966
Punto evaluado:[ 20.76992987 467.28391989  18.25132361]-- Valor de la función en el punto: 55.90864926221509
Punto evaluado:[ 21.51127643 468.4935713   26.83596618]-- Valor de la función en el punto: 56.3079803349592
Punto evaluado:[ 25.10407543 461.50735351  22.57318044]-- Valor de la función en el punto: 53.70658272036371
Punto evaluado:[ 25.12726044 462.02720743  21.21429798]-- Valor de la función en el punto: 53.68813644261122
Punto evaluado:[ 25.1196935  461.58142364  21.76634755]-- Valor de la función en el punto: 53.59729981576747
Punto evaluado:[ 25.11932819 461.55895552  21.73104426]-- Valor de la función en el punto: 53.597307869603654
Punto evaluado:[ 25.11933676 461.55693182  21.73078482]-- Valor de la función en el punto: 53.59733792691605
Punto evaluado:[ 25.1193383  461.5567711   21.73080493]-- Valor de la función en el punto: 53.597340396759385
Punto evaluado:[ 25.11933843 461.55675802  21.73080656]-- Valor de la función en el punto: 53.5973405984328
Punto evaluado:[ 25.11933844 461.55675695  21.73080669]-- Valor de la función en el punto: 53.59734061485335
Punto evaluado:[ 25.11933844 461.55675686  21.7308067 ]-- Valor de la función en el punto: 53.59734061619003

Algoritmos evolutivos

En 2 dimensiones

# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0

# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5

# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)

# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initCycle, creator.Individuo,
                 (toolbox.attr_float, toolbox.attr_float), n=1)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)

# Definir la función de evaluación
def evaluar(individuo):
    return f_griewank2d(individuo[0], individuo[1]),

# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)

# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)

# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
    # Seleccionar a los padres
    padres = toolbox.select(poblacion, len(poblacion))
    # Clonar los padres para aplicar operadores genéticos
    hijos = list(map(toolbox.clone, padres))

    # Aplicar cruce y mutación a los hijos
    for hijo in hijos:
        if random.random() < probabilidad_cruce:
            toolbox.mate(hijo, random.choice(hijos))
            del hijo.fitness.values

        if random.random() < probabilidad_mutacion:
            toolbox.mutate(hijo)
            del hijo.fitness.values

    # Evaluar la aptitud de los nuevos individuos
    individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
    aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
    for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
        individuo.fitness.values = aptitud

    # Reemplazar la población anterior por la nueva generación
    poblacion = hijos

# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])

Mejor solución encontrada: [3.3199442668675464e-06, -2.4718480383282513e-06]
Valor de la función objetivo: 7.042810779012143e-12

c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning:

A class named 'FitnessMin' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.

c:\Users\Hp\anaconda3\envs\tf_gpu\lib\site-packages\deap\creator.py:138: RuntimeWarning:

A class named 'Individuo' has already been created and it will be overwritten. Consider deleting previous creation of that class or rename it.

En 3 dimensiones

# Definir los límites de los parámetros
limite_inferior = -5.0
limite_superior = 5.0

# Definir los parámetros del algoritmo genético
num_generaciones = 50
tamano_poblacion = 100
probabilidad_mutacion = 0.2
probabilidad_cruce = 0.5

# Crear los tipos de objeto para el problema
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individuo", list, fitness=creator.FitnessMin)

# Crear la caja de herramientas del algoritmo genético
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, limite_inferior, limite_superior)
toolbox.register("individuo", tools.initRepeat, creator.Individuo,
                 toolbox.attr_float, n=3)
toolbox.register("poblacion", tools.initRepeat, list, toolbox.individuo)

# Definir la función de evaluación
def evaluar(individuo):
    return f_griewank3d(individuo[0], individuo[1], individuo[2]),

# Definir los operadores genéticos
toolbox.register("evaluate", evaluar)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)

# Crear la población inicial
poblacion = toolbox.poblacion(n=tamano_poblacion)

# Ejecutar el algoritmo genético
for generacion in range(num_generaciones):
    # Seleccionar a los padres
    padres = toolbox.select(poblacion, len(poblacion))
    # Clonar los padres para aplicar operadores genéticos
    hijos = list(map(toolbox.clone, padres))

    # Aplicar cruce y mutación a los hijos
    for hijo in hijos:
        if random.random() < probabilidad_cruce:
            toolbox.mate(hijo, random.choice(hijos))
            del hijo.fitness.values

        if random.random() < probabilidad_mutacion:
            toolbox.mutate(hijo)
            del hijo.fitness.values

    # Evaluar la aptitud de los nuevos individuos
    individuos_nuevos = [individuo for individuo in hijos if not individuo.fitness.valid]
    aptitudes_nuevas = toolbox.map(toolbox.evaluate, individuos_nuevos)
    for individuo, aptitud in zip(individuos_nuevos, aptitudes_nuevas):
        individuo.fitness.values = aptitud

    # Reemplazar la población anterior por la nueva generación
    poblacion = hijos

# Seleccionar al mejor individuo
mejor_individuo = tools.selBest(poblacion, 1)[0]
print("Mejor solución encontrada:", mejor_individuo)
print("Valor de la función objetivo:", mejor_individuo.fitness.values[0])

Mejor solución encontrada: [3.307304896771483e-07, 1.4990202736401472e-06, -6.790513197076734e-07]
Valor de la función objetivo: 6.940004126931854e-13

Optimizacion de particulas

En 2 dimensiones

#Definimos la funcion Griewank para los calculos
def f(x, y):
    return (x ** 2/(4000) + y ** 2/(4000)) - ((np.cos(x/np.sqrt(1))+1) * (np.cos(y/np.sqrt(2))+1))


#Calculamos los datos X Y de la funcion
x, y = np.array(np.meshgrid(np.linspace(-600, 600, 100),
                np.linspace(-600, 600, 100)))
#Evalamos los datos X Y
z = f(x, y)

#Encontramos el minimo Global de la funcion
x_min = x.ravel()[z.argmin()]
y_min = y.ravel()[z.argmin()]
# Se inicializan los Hiperparametros del algoritmo PSO
c1 = c2 = 0.1  # Factores de Acelearacion
w = 0.8  # Coeficiente de Inercia

# Creamos las particulas
n_particles = 35
np.random.seed(100)
X = np.random.rand(2, n_particles) * 500  # Posiscion Inicial de las particulas
V = np.random.randn(2, n_particles) * 200  # velocidad de las particulas

# Inicializacion de datos
pbest = X  # Mejor Posicion de la particula
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()]  # Mejor Posicion del Enjambre
gbest_obj = pbest_obj.min()
#funciond de actualizacion de los puntos del enjambre


def update():
    global V, X, pbest, pbest_obj, gbest, gbest_obj
    # actualizacion de parametros
    r1, r2 = np.random.rand(2)
    # Actualizacion de velocidades
    V = w * V + c1*r1*(pbest - X) + c2*r2*(gbest.reshape(-1, 1)-X)
    X = X + V  # Actualizacion de Posicion de particulas
    obj = f(X[0], X[1])  # Nueva Posicion
    # Se evalua si la posicion nueva 'obj' es mejor a la anterior
    pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]
    pbest_obj = np.array([pbest_obj, obj]).min(axis=0)
    gbest = pbest[:, pbest_obj.argmin()]
    gbest_obj = pbest_obj.min()


# Visualizacion del grafico de Grienwank Y los puntos del enjambre
fig, ax = plt.subplots(figsize=(16, 8))
fig.set_tight_layout(True)
img = ax.imshow(z, extent=[-600, 600, -600, 600],
                origin='lower', cmap='plasma', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([x_min], [y_min], marker='x', markersize=10,
        color="white")  # Minimo Global
contours = ax.contour(x, y, z, 10, colors='black', alpha=0.4)
ax.clabel(contours, inline=True, fontsize=12, fmt="%.0f")
# mejor posision del punto
pbest_plot = ax.scatter(
    pbest[0], pbest[1], marker='o', color='black', alpha=0.5)
p_plot = ax.scatter(X[0], X[1], marker='o', color='green', alpha=0.5)  # Puntos
p_arrow = ax.quiver(X[0], X[1], V[0], V[1], color='green', width=0.005,
                    angles='xy', scale_units='xy', scale=1)  # flecha de los puntos
gbest_plot = plt.scatter([gbest[0]], [gbest[1]], marker='*',
                         s=100, color='red', alpha=0.4)  # Mejor posicion del enjambre
ax.set_xlim([-600, 600])
ax.set_ylim([-600, 600])


def animate(i):
    # Animacion de los puntos en moviento
    title = 'Iteration {:02d}'.format(i)
    # Update params
    update()
    # Set picture
    ax.set_title(title)
    pbest_plot.set_offsets(pbest.T)
    p_plot.set_offsets(X.T)
    p_arrow.set_offsets(X.T)
    p_arrow.set_UVC(V[0], V[1])
    gbest_plot.set_offsets(gbest.reshape(1, -1))
    return ax, pbest_plot, p_plot, p_arrow, gbest_plot


anim = FuncAnimation(fig, animate, frames=list(
    range(1, 50)), interval=500, blit=False, repeat=True)
# Se guarda la animacion en formato gif
anim.save("PSOGRI.gif", dpi=120, writer="imagemagick")
# Informacion genreal del algoritmo luego de optimizar
print("PSO found best solution at f({})={}".format(gbest, gbest_obj))
print("Global optimal at f({})={}".format([x_min,y_min], f(x_min,y_min)))

MovieWriter imagemagick unavailable; using Pillow instead.

PSO found best solution at f([-6.33640917  0.14168459])=-3.977103581529073
Global optimal at f([-6.060606060606119, 18.18181818181813])=-3.776287411893872

Gif Funcion Griewnak

from IPython.display import Image, display

Image(filename='PSOGRI.gif')

<IPython.core.display.Image object>

En 3 dimensiones

lb = [-10, -10, -10]  # límite inferior de las variables
ub = [10, 10, 10]     # límite superior de las variables

x_opt, f_opt = pso(lambda x: f_griewank3d(x[0], x[1], x[2]), lb, ub)

print("El punto óptimo es:", x_opt)

Stopping search: Swarm best objective change less than 1e-08
El punto óptimo es: [-6.27996580e+00  5.09301382e-05  7.34195399e-05]

Optimizacion por evolucion diferencial

En 2 dimensiones

# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones

# Definir el límite inferior y superior para los valores de x e y
lim_inf = -10
lim_sup = 10

# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
    x = random.uniform(lim_inf, lim_sup)
    y = random.uniform(lim_inf, lim_sup)
    pop.append((x, y))

# Realizar la evolución diferencial
for gen in range(ngen):
    for i in range(NP):
        # Seleccionar 3 individuos de forma aleatoria
        a, b, c = random.sample(pop, 3)

        # Seleccionar un número aleatorio para cada variable
        R1, R2 = random.random(), random.random()

        # Realizar la mutación diferencial
        x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
        y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)

        # Realizar la recombinación
        if random.random() < CR:
            x = x_new
            y = y_new
        else:
            x = a[0]
            y = a[1]

        # Evaluar la función en el nuevo individuo
        fitness_new = f_griewank2d(x, y)

        # Reemplazar el individuo anterior si el nuevo individuo es mejor
        if fitness_new < f_griewank2d(pop[i][0], pop[i][1]):
            pop[i] = (x, y)

# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_griewank2d(x[0], x[1]))
print('Valor mínimo de la función: ', f_griewank2d(best_ind[0], best_ind[1]))
print('Variables que producen el valor mínimo: ', best_ind)

Valor mínimo de la función:  0.007396040334114895
Variables que producen el valor mínimo:  (3.140022634806684, 4.438444483718577)

En 3 dimensiones

# Definir los parámetros de la evolución diferencial
F = 0.5 # factor de amplificación
CR = 0.9 # tasa de recombinación
NP = 20 # tamaño de la población
ngen = 100 # número de generaciones

# Definir el límite inferior y superior para los valores de x, y, y z
lim_inf = -10
lim_sup = 10

# Generar la población inicial de forma aleatoria
pop = []
for i in range(NP):
    x = random.uniform(lim_inf, lim_sup)
    y = random.uniform(lim_inf, lim_sup)
    z = random.uniform(lim_inf, lim_sup)
    pop.append((x, y, z))

# Realizar la evolución diferencial
for gen in range(ngen):
    for i in range(NP):
        # Seleccionar 3 individuos de forma aleatoria
        a, b, c = random.sample(pop, 3)

        # Seleccionar un número aleatorio para cada variable
        R1, R2, R3 = random.random(), random.random(), random.random()

        # Realizar la mutación diferencial
        x_new = np.clip(a[0] + F*(b[0] - c[0]), lim_inf, lim_sup)
        y_new = np.clip(a[1] + F*(b[1] - c[1]), lim_inf, lim_sup)
        z_new = np.clip(a[2] + F*(b[2] - c[2]), lim_inf, lim_sup)

        # Realizar la recombinación
        if random.random() < CR:
            x = x_new
            y = y_new
            z = z_new
        else:
            x = a[0]
            y = a[1]
            z = a[2]

        # Evaluar la función en el nuevo individuo
        fitness_new = f_griewank3d(x, y, z)

        # Reemplazar el individuo anterior si el nuevo individuo es mejor
        if fitness_new < f_griewank3d(pop[i][0], pop[i][1], pop[i][2]):
            pop[i] = (x, y, z)

# Imprimir el mejor individuo encontrado
best_ind = min(pop, key=lambda x: f_griewank3d(x[0], x[1], x[2]))
print('Valor mínimo de la función: ', f_griewank3d(best_ind[0], best_ind[1], best_ind[2]))
print('Variables que producen el valor mínimo: ', best_ind)

Valor mínimo de la función:  0.056348699247850464
Variables que producen el valor mínimo:  (3.062545509753079, 0.38435240630277434, -5.639937371581347)

Parte 2.1: Problema Viajero Colombian Version.

#Importamos las librerias Necesarias para el desarollo Utlizando Algoritmos Geneticos
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Point
import pyproj 
from pygad import gann
from matplotlib import animation, rc
import pandas as pd
import random
import itertools
import imageio
import os

Clase AntColony para el desarrollo del Ejercicio.

import numpy as np
import matplotlib.pyplot as plt
import time

import warnings

warnings.filterwarnings("ignore")


class AntColonyOptimizer:
    def __init__(self, ants, evaporation_rate, intensification, alpha=1.0, beta=0.0, beta_evaporation_rate=0,
                 choose_best=.1):
        """
        Ant colony optimizer.  Traverses a graph and finds either the max or min distance between nodes.
        :param ants: number of ants to traverse the graph
        :param evaporation_rate: rate at which pheromone evaporates
        :param intensification: constant added to the best path
        :param alpha: weighting of pheromone
        :param beta: weighting of heuristic (1/distance)
        :param beta_evaporation_rate: rate at which beta decays (optional)
        :param choose_best: probability to choose the best route
        """
        # Parameters
        self.ants = ants
        self.evaporation_rate = evaporation_rate
        self.pheromone_intensification = intensification
        self.heuristic_alpha = alpha
        self.heuristic_beta = beta
        self.beta_evaporation_rate = beta_evaporation_rate
        self.choose_best = choose_best

        # Internal representations
        self.pheromone_matrix = None
        self.heuristic_matrix = None
        self.probability_matrix = None

        self.map = None
        self.set_of_available_nodes = None

        # Internal stats
        self.best_series = []
        self.best = None
        self.fitted = False
        self.best_path = None
        self.fit_time = None

        # Plotting values
        self.stopped_early = False

    def __str__(self):
        string = "Ant Colony Optimizer"
        string += "\n--------------------"
        string += "\nDesigned to optimize either the minimum or maximum distance between nodes in a square matrix that behaves like a distance matrix."
        string += "\n--------------------"
        string += "\nNumber of ants:\t\t\t\t{}".format(self.ants)
        string += "\nEvaporation rate:\t\t\t{}".format(self.evaporation_rate)
        string += "\nIntensification factor:\t\t{}".format(
            self.pheromone_intensification)
        string += "\nAlpha Heuristic:\t\t\t{}".format(self.heuristic_alpha)
        string += "\nBeta Heuristic:\t\t\t\t{}".format(self.heuristic_beta)
        string += "\nBeta Evaporation Rate:\t\t{}".format(
            self.beta_evaporation_rate)
        string += "\nChoose Best Percentage:\t\t{}".format(self.choose_best)
        string += "\n--------------------"
        string += "\nUSAGE:"
        string += "\nNumber of ants influences how many paths are explored each iteration."
        string += "\nThe alpha and beta heuristics affect how much influence the pheromones or the distance heuristic weigh an ants' decisions."
        string += "\nBeta evaporation reduces the influence of the heuristic over time."
        string += "\nChoose best is a percentage of how often an ant will choose the best route over probabilistically choosing a route based on pheromones."
        string += "\n--------------------"
        if self.fitted:
            string += "\n\nThis optimizer has been fitted."
        else:
            string += "\n\nThis optimizer has NOT been fitted."
        return string

    def _initialize(self):
        """
        Initializes the model by creating the various matrices and generating the list of available nodes
        """
        assert self.map.shape[0] == self.map.shape[1], "Map is not a distance matrix!"
        num_nodes = self.map.shape[0]
        self.pheromone_matrix = np.ones((num_nodes, num_nodes))
        # Remove the diagonal since there is no pheromone from node i to itself
        self.pheromone_matrix[np.eye(num_nodes) == 1] = 0
        self.heuristic_matrix = 1 / self.map
        self.probability_matrix = (self.pheromone_matrix ** self.heuristic_alpha) * (
            self.heuristic_matrix ** self.heuristic_beta)  # element by element multiplcation
        self.set_of_available_nodes = list(range(num_nodes))

    def _reinstate_nodes(self):
        """
        Resets available nodes to all nodes for the next iteration
        """
        self.set_of_available_nodes = list(range(self.map.shape[0]))

    def _update_probabilities(self):
        """
        After evaporation and intensification, the probability matrix needs to be updated.  This function
        does that.
        """
        self.probability_matrix = (self.pheromone_matrix ** self.heuristic_alpha) * (
            self.heuristic_matrix ** self.heuristic_beta)

    def _choose_next_node(self, from_node):
        """
        Chooses the next node based on probabilities.  If p < p_choose_best, then the best path is chosen, otherwise
        it is selected from a probability distribution weighted by the pheromone.
        :param from_node: the node the ant is coming from
        :return: index of the node the ant is going to
        """
        numerator = self.probability_matrix[from_node,
                                            self.set_of_available_nodes]
        if np.random.random() < self.choose_best:
            next_node = np.argmax(numerator)
        else:
            denominator = np.sum(numerator)
            probabilities = numerator / denominator
            next_node = np.random.choice(
                range(len(probabilities)), p=probabilities)
        return next_node

    def _remove_node(self, node):
        self.set_of_available_nodes.remove(node)

    def _evaluate(self, paths, mode):
        """
        Evaluates the solutions of the ants by adding up the distances between nodes.
        :param paths: solutions from the ants
        :param mode: max or min
        :return: x and y coordinates of the best path as a tuple, the best path, and the best score
        """
        scores = np.zeros(len(paths))
        coordinates_i = []
        coordinates_j = []
        for index, path in enumerate(paths):
            score = 0
            coords_i = []
            coords_j = []
            for i in range(len(path) - 1):
                coords_i.append(path[i])
                coords_j.append(path[i + 1])
                score += self.map[path[i], path[i + 1]]
            scores[index] = score
            coordinates_i.append(coords_i)
            coordinates_j.append(coords_j)
        if mode == 'min':
            best = np.argmin(scores)
        elif mode == 'max':
            best = np.argmax(scores)
        return (coordinates_i[best], coordinates_j[best]), paths[best], scores[best]

    def _evaporation(self):
        """
        Evaporate some pheromone as the inverse of the evaporation rate.  Also evaporates beta if desired.
        """
        self.pheromone_matrix *= (1 - self.evaporation_rate)
        self.heuristic_beta *= (1 - self.beta_evaporation_rate)

    def _intensify(self, best_coords):
        """
        Increases the pheromone by some scalar for the best route.
        :param best_coords: x and y (i and j) coordinates of the best route
        """
        i = best_coords[0]
        j = best_coords[1]
        self.pheromone_matrix[i, j] += self.pheromone_intensification

    def fit(self, map_matrix, iterations=100, mode='min', early_stopping_count=20, verbose=True):
        """
        Fits the ACO to a specific map.  This was designed with the Traveling Salesman problem in mind.
        :param map_matrix: Distance matrix or some other matrix with similar properties
        :param iterations: number of iterations
        :param mode: whether to get the minimum path or maximum path
        :param early_stopping_count: how many iterations of the same score to make the algorithm stop early
        :return: the best score
        """
        if verbose:
            print("Beginning ACO Optimization with {} iterations...".format(iterations))
        self.map = map_matrix
        start = time.time()
        self._initialize()
        num_equal = 0

        for i in range(iterations):
            start_iter = time.time()
            paths = []
            path = []

            for ant in range(self.ants):
                current_node = self.set_of_available_nodes[np.random.randint(
                    0, len(self.set_of_available_nodes))]
                start_node = current_node
                while True:
                    path.append(current_node)
                    self._remove_node(current_node)
                    if len(self.set_of_available_nodes) != 0:
                        current_node_index = self._choose_next_node(
                            current_node)
                        current_node = self.set_of_available_nodes[current_node_index]
                    else:
                        break

                path.append(start_node)  # go back to start
                self._reinstate_nodes()
                paths.append(path)
                path = []

            best_path_coords, best_path, best_score = self._evaluate(
                paths, mode)

            if i == 0:
                best_score_so_far = best_score
            else:
                if mode == 'min':
                    if best_score < best_score_so_far:
                        best_score_so_far = best_score
                        self.best_path = best_path
                elif mode == 'max':
                    if best_score > best_score_so_far:
                        best_score_so_far = best_score
                        self.best_path = best_path

            if best_score == best_score_so_far:
                num_equal += 1
            else:
                num_equal = 0

            self.best_series.append(best_score)
            self._evaporation()
            self._intensify(best_path_coords)
            self._update_probabilities()

            if verbose:
                print("Best score at iteration {}: {}; overall: {} ({}s)"
                      "".format(i, round(best_score, 2), round(best_score_so_far, 2),
                                round(time.time() - start_iter)))

            if best_score == best_score_so_far and num_equal == early_stopping_count:
                self.stopped_early = True
                print("Stopping early due to {} iterations of the same score.".format(
                    early_stopping_count))
                break

        self.fit_time = round(time.time() - start)
        self.fitted = True

        if mode == 'min':
            self.best = self.best_series[np.argmin(self.best_series)]
            if verbose:
                print(
                    "ACO fitted.  Runtime: {} minutes.  Best score: {}".format(self.fit_time // 60, self.best))
            return self.best
        elif mode == 'max':
            self.best = self.best_series[np.argmax(self.best_series)]
            if verbose:
                print(
                    "ACO fitted.  Runtime: {} minutes.  Best score: {}".format(self.fit_time // 60, self.best))
            return self.best
        else:
            raise ValueError("Invalid mode!  Choose 'min' or 'max'.")

    def plot(self):
        """
        Plots the score over time after the model has been fitted.
        :return: None if the model isn't fitted yet
        """
        if not self.fitted:
            print("Ant Colony Optimizer not fitted!  There exists nothing to plot.")
            return None
        else:
            fig, ax = plt.subplots(figsize=(20, 15))
            ax.plot(self.best_series, label="Best Run")
            ax.set_xlabel("Iteration")
            ax.set_ylabel("Performance")
            ax.text(.8, .6,
                    'Ants: {}\nEvap Rate: {}\nIntensify: {}\nAlpha: {}\nBeta: {}\nBeta Evap: {}\nChoose Best: {}\n\nFit Time: {}m{}'.format(
                        self.ants, self.evaporation_rate, self.pheromone_intensification, self.heuristic_alpha,
                        self.heuristic_beta, self.beta_evaporation_rate, self.choose_best, self.fit_time // 60,
                        ["\nStopped Early!" if self.stopped_early else ""][0]),
                    bbox={'facecolor': 'gray', 'alpha': 0.8, 'pad': 10}, transform=ax.transAxes)
            ax.legend()
            plt.title("Ant Colony Optimization Results (best: {})".format(
                np.round(self.best, 2)))
            plt.show()

Instanciamos las Ciudades y las coordenadas

# Coordenades X Y de las ciudades
ciudades = ["Palmira", "Pasto", "Tuluá", "Bogota", "Pereira", "Armenia", "Manizales", "Valledupar",
"Montería", "Soledad", "Cartagena", "Barranquilla", "Medellín", "Bucaramanga", "Cúcuta"]

idx_ciudades=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14]

#Las coordenadas corresponden al orden en como esta la lista de ciudades
coords=[[3.53944, -76.30361], [1.21361, -77.28111], [4.08466, -76.19536], [4.60971, -74.08175], [4.81333, -75.69611], [4.53389, -75.68111], [5.06889, -75.51738], [10.46314, -73.25322], 
[8.74798, -75.88143], [10.91843, -74.76459], [10.39972, -75.51444], [10.96854, -74.78132], [6.25184, -75.56359], [7.12539, -73.1198], [7.89391, -72.50782]]

Y = [coords[i][0] for i in range(len(coords))]
X = [coords[i][1] for i in range(len(coords))]
print(X)
print(Y)

[-76.30361, -77.28111, -76.19536, -74.08175, -75.69611, -75.68111, -75.51738, -73.25322, -75.88143, -74.76459, -75.51444, -74.78132, -75.56359, -73.1198, -72.50782]
[3.53944, 1.21361, 4.08466, 4.60971, 4.81333, 4.53389, 5.06889, 10.46314, 8.74798, 10.91843, 10.39972, 10.96854, 6.25184, 7.12539, 7.89391]

Se Realiza una conversion de loas coordenadas Anteriores para que coincidan con el formato del archvo .shp del mapa de Colombia

#TRANSFORMACION DE LAS COORDENAS
# Definir el sistema de proyección de origen (WGS84)
source_crs = 'EPSG:4326'

# Definir el sistema de proyección de destino 
target_crs = 'EPSG:21897'
# Crear un objeto de transformación de proyección
transformer = pyproj.Transformer.from_crs(
    source_crs, target_crs, always_xy=True)
X_T, Y_T = transformer.transform(X, Y)
print(X_T)
print(Y_T)

# Crear una lista de puntos Shapely a partir de las coordenadas X e Y
geometry = [Point(x, y) for x, y in zip(X_T, Y_T)]
# Crear un GeoDataFrame con las ciudades y la geometría
df = gpd.GeoDataFrame(
    {'Ciudad': ciudades, 'geometry': geometry})

#Cambia el formato de POINT a FLOAT puesto que mas adelanto realizaremos calculos con estos
def punto_a_entero(punto):
    x, y = punto.x, punto.y
    return float(x), float(y)
#Agregamos la columna Punto_float al dataframe
df['punto_float'] = df['geometry'].apply(lambda punto: punto_a_entero(punto))
df

[752605.9915458227, 643275.0465845948, 764787.6710086607, 999529.3630757327, 820426.0323943375, 822021.5505042156, 840318.4461557847, 1090247.024864841, 801474.1002934285, 924887.5034194184, 842633.5437126327, 923071.5071483594, 835526.5527440282, 1105793.363446178, 1173108.296118954]
[883431.2205430069, 626167.8831370582, 943734.5434153986, 1001493.9629269837, 1024217.0013714742, 993300.2959815342, 1052441.7478454192, 1648963.4840482422, 1459613.1832835586, 1699295.6026588602, 1642191.7124144642, 1704843.1924954015, 1183313.114715379, 1279803.6634820616, 1365014.5786278923]

	Ciudad	geometry	punto_float
0	Palmira	POINT (752605.992 883431.221)	(752605.9915458227, 883431.2205430069)
1	Pasto	POINT (643275.047 626167.883)	(643275.0465845948, 626167.8831370582)
2	Tuluá	POINT (764787.671 943734.543)	(764787.6710086607, 943734.5434153986)
3	Bogota	POINT (999529.363 1001493.963)	(999529.3630757327, 1001493.9629269837)
4	Pereira	POINT (820426.032 1024217.001)	(820426.0323943375, 1024217.0013714742)
5	Armenia	POINT (822021.551 993300.296)	(822021.5505042156, 993300.2959815342)
6	Manizales	POINT (840318.446 1052441.748)	(840318.4461557847, 1052441.7478454192)
7	Valledupar	POINT (1090247.025 1648963.484)	(1090247.024864841, 1648963.4840482422)
8	Montería	POINT (801474.100 1459613.183)	(801474.1002934285, 1459613.1832835586)
9	Soledad	POINT (924887.503 1699295.603)	(924887.5034194184, 1699295.6026588602)
10	Cartagena	POINT (842633.544 1642191.712)	(842633.5437126327, 1642191.7124144642)
11	Barranquilla	POINT (923071.507 1704843.192)	(923071.5071483594, 1704843.1924954015)
12	Medellín	POINT (835526.553 1183313.115)	(835526.5527440282, 1183313.114715379)
13	Bucaramanga	POINT (1105793.363 1279803.663)	(1105793.363446178, 1279803.6634820616)
14	Cúcuta	POINT (1173108.296 1365014.579)	(1173108.296118954, 1365014.5786278923)

# Cargar el shapefile de Colombia
mapa = gpd.read_file('Colombia Shape\depto.shp')
#print(mapa.crs) #Formato CRS del archivo .shp
fig, ax = plt.subplots(figsize=(14, 14))
mapa.plot(ax=ax, color='lime', edgecolor='black')
df.plot(ax=ax, markersize=50, color='red')
for i, txt in enumerate(df['Ciudad']):
    ax.annotate(txt, (df['geometry'][i].x,df['geometry'][i].y), fontsize=12)

Matriz De Costos.

Matrices Requeridas:

MATRIZ DE DISTANCIAS ENTRE LAS CIUDADES:MDC

MATRIZ DE # DE PEJAES ENTRE LAS CIUDADES:MNP

MATRIZ DE TIEMPO ESTIMADO ENTRE CIUDADES:MTC

Para realizar estos calculos tambien hay que tener en cuenta la informacion del Automovil que se empleara para este problema el cual es un Chevrolet Camaro SS Este Camaro usa Gasolina como combustible (Inserte Imagen de Camaro)

De este Auto tenemos lo siguiente:

Su rendimeinto en KM por Galon es: 29.15Km/Galon.

El precio de la gasolina a la fecha de realizacion de este trabajo (Marzo 15 del 2023) es de $9.664 Pesos Colombianos por Galon.

EL sueldo por Hora de nuestro chofer a la fecha de realizacion de este trabajo (Marzo 15 del 2023) sera de $1368,39 COP

Nuestro automovil en Colombia pagaria la tarifa de categoria 1. Este costo sera de Aproximadamente $8000 COP

Matriz de Distancias Entre las ciudades: MDC

Para esto Usaremos la Formula de Haversine que sirve para calcular la distancia en km entre 2 puntos de la tierra el cual nos dara un estimado de las diatancias que requerimos para la solucion de este problema.

#FUNCION DE HAVERSINE
from math import radians, sin, cos, sqrt, atan2
def haversine(coord1, coord2):
    # Radio de la Tierra en km
    R = 6371.0

    # Convertir coordenadas de grados a radianes
    lat1, lon1 = map(radians, coord1)
    lat2, lon2 = map(radians, coord2)

    # Diferencia de latitud y longitud
    dlat = lat2 - lat1
    dlon = lon2 - lon1

    # Fórmula de Haversine
    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
    c = 2 * atan2(sqrt(a), sqrt(1 - a))

    # Distancia entre las dos coordenadas
    distance = R * c

    return distance

# Creamos una matriz vacía de distancias
distances = [[0.0] * len(coords) for i in range(len(coords))]

# Calcular la distancia entre cada par de ciudades
for i in range(len(coords)):
    for j in range(i + 1, len(coords)):
        # Calcular la distancia entre las coordenadas
        dist = haversine(coords[i], coords[j])

        # Almacenar la distancia en la matriz y los redondeamos a 2 decimales
        distances[i][j] = round(dist,2)
        distances[j][i] = round(dist,2)

# Imprimir la matriz de distancias
#for i in range(len(ciudades)):
 #   print(ciudades[i], distances[i])
#convertimos la lista en una matriz 15X15
MDC = np.array(distances).reshape(15,15) 
MDC

array([[   0.  ,  280.49,   61.8 ,  273.66,  156.85,  130.36,  191.11,
         840.18,  581.04,  837.86,  767.78,  842.95,  312.55,  532.17,
         640.89],
       [ 280.49,    0.  ,  341.26,  518.46,  437.24,  409.72,  471.28,
        1120.66,  851.99, 1114.34, 1039.94, 1119.27,  591.73,  803.06,
         911.8 ],
       [  61.8 ,  341.26,    0.  ,  241.51,   98.12,   75.81,  132.76,
         779.9 ,  519.7 ,  776.06,  706.2 ,  781.14,  250.93,  479.72,
         587.92],
       [ 273.66,  518.46,  241.51,    0.  ,  180.33,  177.48,  167.06,
         657.24,  501.23,  705.52,  662.89,  711.25,  245.45,  299.28,
         404.5 ],
       [ 156.85,  437.24,   98.12,  180.33,    0.  ,   31.12,   34.63,
         683.44,  437.99,  686.56,  621.5 ,  691.8 ,  160.63,  383.75,
         491.38],
       [ 130.36,  409.72,   75.81,  177.48,   31.12,    0.  ,   62.19,
         711.52,  469.11,  717.07,  652.51,  722.33,  191.47,  404.1 ,
         512.45],
       [ 191.11,  471.28,  132.76,  167.06,   34.63,   62.19,    0.  ,
         649.58,  411.06,  655.7 ,  592.76,  660.99,  131.64,  350.08,
         457.4 ],
       [ 840.18, 1120.66,  779.9 ,  657.24,  683.44,  711.52,  649.58,
           0.  ,  345.53,  172.72,  247.38,  176.16,  532.78,  371.43,
         297.17],
       [ 581.04,  851.99,  519.7 ,  501.23,  437.99,  469.11,  411.06,
         345.53,    0.  ,  270.59,  188.02,  274.75,  279.76,  353.62,
         383.13],
       [ 837.86, 1114.34,  776.06,  705.52,  686.56,  717.07,  655.7 ,
         172.72,  270.59,    0.  ,  100.2 ,    5.86,  526.28,  458.8 ,
         417.59],
       [ 767.78, 1039.94,  706.2 ,  662.89,  621.5 ,  652.51,  592.76,
         247.38,  188.02,  100.2 ,    0.  ,  102.07,  461.25,  449.22,
         431.93],
       [ 842.95, 1119.27,  781.14,  711.25,  691.8 ,  722.33,  660.99,
         176.16,  274.75,    5.86,  102.07,    0.  ,  531.47,  464.64,
         423.16],
       [ 312.55,  591.73,  250.93,  245.45,  160.63,  191.47,  131.64,
         532.78,  279.76,  526.28,  461.25,  531.47,    0.  ,  286.83,
         383.45],
       [ 532.17,  803.06,  479.72,  299.28,  383.75,  404.1 ,  350.08,
         371.43,  353.62,  458.8 ,  449.22,  464.64,  286.83,    0.  ,
         108.88],
       [ 640.89,  911.8 ,  587.92,  404.5 ,  491.38,  512.45,  457.4 ,
         297.17,  383.13,  417.59,  431.93,  423.16,  383.45,  108.88,
           0.  ]])

Ahora Procedemos a las creacion de las matrices MNP y MTC las que corresponden al numero de peajes entre ciudades y el tiempo para viajar de una ciudad a otra

#procedemos a convertir los archivos .xlsx a matrices 
# Cargar el archivo de Excel en un DataFrame
Peajes = pd.read_excel('PeajesSimples.xlsx',header=None)

# Convertir el DataFrame en un array de NumPy
MNP = Peajes.to_numpy()

Tiempos = pd.read_excel('TiemposSimples.xlsx', header=None) #en Total de Minutos

# Convertir el DataFrame en un array de NumPy
MTC = Tiempos.to_numpy()

print(MNP.shape)
print(MTC.shape)

(15, 15)
(15, 15)

Con estas 3 matrices ya calculadas damos luz para la matriz de costos totales de movilizacion entre ciudades la denotaremos MCTC

#Datos usados para los caclculos
Vh = 1368.39 # Salrio de una hora de sueldo del chofer
Cpc1 = 8000 # Costo Aproximado de los peajes (Esto puede variar mucho o poco entre peajes para la practicidad del ejercicio se fijo).
Pg = 9.664 # PRecio del Galon de Combustible

#Realizamos alguno Ajustes en las matrices 
Txh = MTC/60 # Convertimos los minutos en horas
Dxg = MDC/29.15 # Con esta operacion conseguimos el numero de kilometros por galon 

#FINALMENTE Calculamos la matriz MCTC

MCTC = Vh*Txh + MNP*8000 + Pg*Dxg 

print(MCTC.shape)

(15, 15)

Problema del Viajero Colombian Version: Colonias De Hormigas

problem = MCTC # matriz de costos totales
Hormigas=[10,20,80,100,120] # nuemro de Hormigas a evaluar 
MejorRuta=[] #Guardamos la ruta
for i in range(len(Hormigas)):
    
    optimizer = AntColonyOptimizer(ants=Hormigas[i], evaporation_rate=.1, intensification=2, alpha=1, beta=1,
                                beta_evaporation_rate=0.5, choose_best=.1)
    best = optimizer.fit(problem, 100) # Iniciamos el optimizador con 100 iteraciones
    optimizer.plot()
    MejorRuta.append(optimizer.best_path) #Agregamos la mejor ruta por cada numero de Hormigas

Beginning ACO Optimization with 100 iterations...
Best score at iteration 0: 958526.71; overall: 958526.71 (0s)
Best score at iteration 1: 833377.94; overall: 833377.94 (0s)
Best score at iteration 2: 998090.98; overall: 833377.94 (0s)
Best score at iteration 3: 1094889.29; overall: 833377.94 (0s)
Best score at iteration 4: 1157027.73; overall: 833377.94 (0s)
Best score at iteration 5: 911612.27; overall: 833377.94 (0s)
Best score at iteration 6: 990091.12; overall: 833377.94 (0s)
Best score at iteration 7: 1063159.76; overall: 833377.94 (0s)
Best score at iteration 8: 878125.67; overall: 833377.94 (0s)
Best score at iteration 9: 1105489.97; overall: 833377.94 (0s)
Best score at iteration 10: 946840.45; overall: 833377.94 (0s)
Best score at iteration 11: 850775.99; overall: 833377.94 (0s)
Best score at iteration 12: 900991.92; overall: 833377.94 (0s)
Best score at iteration 13: 1016693.63; overall: 833377.94 (0s)
Best score at iteration 14: 824877.43; overall: 824877.43 (0s)
Best score at iteration 15: 962244.56; overall: 824877.43 (0s)
Best score at iteration 16: 865798.54; overall: 824877.43 (0s)
Best score at iteration 17: 817484.4; overall: 817484.4 (0s)
Best score at iteration 18: 755048.45; overall: 755048.45 (0s)
Best score at iteration 19: 880523.13; overall: 755048.45 (0s)
Best score at iteration 20: 867293.41; overall: 755048.45 (0s)
Best score at iteration 21: 933027.84; overall: 755048.45 (0s)
Best score at iteration 22: 874751.98; overall: 755048.45 (0s)
Best score at iteration 23: 841491.38; overall: 755048.45 (0s)
Best score at iteration 24: 783664.25; overall: 755048.45 (0s)
Best score at iteration 25: 860620.6; overall: 755048.45 (0s)
Best score at iteration 26: 736701.69; overall: 736701.69 (0s)
Best score at iteration 27: 736701.69; overall: 736701.69 (0s)
Best score at iteration 28: 782859.3; overall: 736701.69 (0s)
Best score at iteration 29: 898906.44; overall: 736701.69 (0s)
Best score at iteration 30: 778318.48; overall: 736701.69 (0s)
Best score at iteration 31: 773751.46; overall: 736701.69 (0s)
Best score at iteration 32: 747642.54; overall: 736701.69 (0s)
Best score at iteration 33: 736701.69; overall: 736701.69 (0s)
Best score at iteration 34: 736701.69; overall: 736701.69 (0s)
Best score at iteration 35: 736701.69; overall: 736701.69 (0s)
Best score at iteration 36: 736701.69; overall: 736701.69 (0s)
Best score at iteration 37: 736701.69; overall: 736701.69 (0s)
Best score at iteration 38: 736701.69; overall: 736701.69 (0s)
Best score at iteration 39: 736701.69; overall: 736701.69 (0s)
Best score at iteration 40: 736701.69; overall: 736701.69 (0s)
Best score at iteration 41: 736701.69; overall: 736701.69 (0s)
Best score at iteration 42: 736701.69; overall: 736701.69 (0s)
Best score at iteration 43: 736701.69; overall: 736701.69 (0s)
Best score at iteration 44: 736701.69; overall: 736701.69 (0s)
Best score at iteration 45: 736701.69; overall: 736701.69 (0s)
Best score at iteration 46: 736701.69; overall: 736701.69 (0s)
Best score at iteration 47: 736701.69; overall: 736701.69 (0s)
Best score at iteration 48: 736701.69; overall: 736701.69 (0s)
Best score at iteration 49: 729205.86; overall: 729205.86 (0s)
Best score at iteration 50: 729205.86; overall: 729205.86 (0s)
Best score at iteration 51: 729205.86; overall: 729205.86 (0s)
Best score at iteration 52: 729205.86; overall: 729205.86 (0s)
Best score at iteration 53: 729205.86; overall: 729205.86 (0s)
Best score at iteration 54: 729205.86; overall: 729205.86 (0s)
Best score at iteration 55: 729205.86; overall: 729205.86 (0s)
Best score at iteration 56: 729205.86; overall: 729205.86 (0s)
Best score at iteration 57: 729205.86; overall: 729205.86 (0s)
Best score at iteration 58: 729205.86; overall: 729205.86 (0s)
Best score at iteration 59: 729205.86; overall: 729205.86 (0s)
Best score at iteration 60: 729205.86; overall: 729205.86 (0s)
Best score at iteration 61: 729205.86; overall: 729205.86 (0s)
Best score at iteration 62: 729205.86; overall: 729205.86 (0s)
Best score at iteration 63: 729205.86; overall: 729205.86 (0s)
Best score at iteration 64: 729205.86; overall: 729205.86 (0s)
Best score at iteration 65: 729205.86; overall: 729205.86 (0s)
Best score at iteration 66: 729205.86; overall: 729205.86 (0s)
Best score at iteration 67: 729205.86; overall: 729205.86 (0s)
Best score at iteration 68: 729205.86; overall: 729205.86 (0s)
Best score at iteration 69: 729205.86; overall: 729205.86 (0s)
Best score at iteration 70: 729205.86; overall: 729205.86 (0s)
Best score at iteration 71: 729205.86; overall: 729205.86 (0s)
Best score at iteration 72: 729205.86; overall: 729205.86 (0s)
Best score at iteration 73: 729205.86; overall: 729205.86 (0s)
Best score at iteration 74: 729205.86; overall: 729205.86 (0s)
Best score at iteration 75: 729205.86; overall: 729205.86 (0s)
Best score at iteration 76: 729205.86; overall: 729205.86 (0s)
Best score at iteration 77: 729205.86; overall: 729205.86 (0s)
Best score at iteration 78: 729205.86; overall: 729205.86 (0s)
Best score at iteration 79: 729205.86; overall: 729205.86 (0s)
Best score at iteration 80: 729205.86; overall: 729205.86 (0s)
Best score at iteration 81: 729205.86; overall: 729205.86 (0s)
Best score at iteration 82: 729205.86; overall: 729205.86 (0s)
Best score at iteration 83: 729205.86; overall: 729205.86 (0s)
Best score at iteration 84: 729205.86; overall: 729205.86 (0s)
Best score at iteration 85: 729205.86; overall: 729205.86 (0s)
Best score at iteration 86: 729205.86; overall: 729205.86 (0s)
Best score at iteration 87: 729205.86; overall: 729205.86 (0s)
Best score at iteration 88: 729205.86; overall: 729205.86 (0s)
Best score at iteration 89: 729205.86; overall: 729205.86 (0s)
Best score at iteration 90: 729205.86; overall: 729205.86 (0s)
Best score at iteration 91: 729205.86; overall: 729205.86 (0s)
Best score at iteration 92: 729205.86; overall: 729205.86 (0s)
Best score at iteration 93: 729205.86; overall: 729205.86 (0s)
Best score at iteration 94: 729205.86; overall: 729205.86 (0s)
Best score at iteration 95: 729205.86; overall: 729205.86 (0s)
Best score at iteration 96: 729205.86; overall: 729205.86 (0s)
Best score at iteration 97: 729205.86; overall: 729205.86 (0s)
Best score at iteration 98: 729205.86; overall: 729205.86 (0s)
Best score at iteration 99: 729205.86; overall: 729205.86 (0s)
ACO fitted.  Runtime: 0 minutes.  Best score: 729205.8646025728

Beginning ACO Optimization with 100 iterations...
Best score at iteration 0: 676258.8; overall: 676258.8 (0s)
Best score at iteration 1: 1115739.0; overall: 676258.8 (0s)
Best score at iteration 2: 848845.0; overall: 676258.8 (0s)
Best score at iteration 3: 980641.06; overall: 676258.8 (0s)
Best score at iteration 4: 929794.16; overall: 676258.8 (0s)
Best score at iteration 5: 815519.3; overall: 676258.8 (0s)
Best score at iteration 6: 937332.46; overall: 676258.8 (0s)
Best score at iteration 7: 761930.4; overall: 676258.8 (0s)
Best score at iteration 8: 833334.04; overall: 676258.8 (0s)
Best score at iteration 9: 880173.73; overall: 676258.8 (0s)
Best score at iteration 10: 737772.55; overall: 676258.8 (0s)
Best score at iteration 11: 700355.48; overall: 676258.8 (0s)
Best score at iteration 12: 806433.02; overall: 676258.8 (0s)
Best score at iteration 13: 800018.72; overall: 676258.8 (0s)
Best score at iteration 14: 853380.28; overall: 676258.8 (0s)
Best score at iteration 15: 737772.55; overall: 676258.8 (0s)
Best score at iteration 16: 767859.24; overall: 676258.8 (0s)
Best score at iteration 17: 684849.0; overall: 676258.8 (0s)
Best score at iteration 18: 720133.72; overall: 676258.8 (0s)
Best score at iteration 19: 701458.38; overall: 676258.8 (0s)
Best score at iteration 20: 684418.67; overall: 676258.8 (0s)
Best score at iteration 21: 737772.55; overall: 676258.8 (0s)
Best score at iteration 22: 709301.21; overall: 676258.8 (0s)
Best score at iteration 23: 737772.55; overall: 676258.8 (0s)
Best score at iteration 24: 735301.04; overall: 676258.8 (0s)
Best score at iteration 25: 700355.48; overall: 676258.8 (0s)
Best score at iteration 26: 685076.97; overall: 676258.8 (0s)
Best score at iteration 27: 685076.97; overall: 676258.8 (0s)
Best score at iteration 28: 708332.91; overall: 676258.8 (0s)
Best score at iteration 29: 735883.45; overall: 676258.8 (0s)
Best score at iteration 30: 684418.67; overall: 676258.8 (0s)
Best score at iteration 31: 685076.97; overall: 676258.8 (0s)
Best score at iteration 32: 676258.8; overall: 676258.8 (0s)
Best score at iteration 33: 737772.55; overall: 676258.8 (0s)
Best score at iteration 34: 684418.67; overall: 676258.8 (0s)
Best score at iteration 35: 668281.37; overall: 668281.37 (0s)
Best score at iteration 36: 668281.37; overall: 668281.37 (0s)
Best score at iteration 37: 684418.67; overall: 668281.37 (0s)
Best score at iteration 38: 684418.67; overall: 668281.37 (0s)
Best score at iteration 39: 676258.8; overall: 668281.37 (0s)
Best score at iteration 40: 668281.37; overall: 668281.37 (0s)
Best score at iteration 41: 668281.37; overall: 668281.37 (0s)
Best score at iteration 42: 668281.37; overall: 668281.37 (0s)
Best score at iteration 43: 668281.37; overall: 668281.37 (0s)
Best score at iteration 44: 668281.37; overall: 668281.37 (0s)
Best score at iteration 45: 668281.37; overall: 668281.37 (0s)
Best score at iteration 46: 668281.37; overall: 668281.37 (0s)
Best score at iteration 47: 668281.37; overall: 668281.37 (0s)
Best score at iteration 48: 668281.37; overall: 668281.37 (0s)
Best score at iteration 49: 668281.37; overall: 668281.37 (0s)
Best score at iteration 50: 668281.37; overall: 668281.37 (0s)
Best score at iteration 51: 668281.37; overall: 668281.37 (0s)
Best score at iteration 52: 668281.37; overall: 668281.37 (0s)
Best score at iteration 53: 668281.37; overall: 668281.37 (0s)
Best score at iteration 54: 668281.37; overall: 668281.37 (0s)
Best score at iteration 55: 668281.37; overall: 668281.37 (0s)
Best score at iteration 56: 668281.37; overall: 668281.37 (0s)
Best score at iteration 57: 668281.37; overall: 668281.37 (0s)
Best score at iteration 58: 668281.37; overall: 668281.37 (0s)
Best score at iteration 59: 668281.37; overall: 668281.37 (0s)
Best score at iteration 60: 668281.37; overall: 668281.37 (0s)
Best score at iteration 61: 668281.37; overall: 668281.37 (0s)
Best score at iteration 62: 668281.37; overall: 668281.37 (0s)
Stopping early due to 20 iterations of the same score.
ACO fitted.  Runtime: 0 minutes.  Best score: 668281.3692114921

Beginning ACO Optimization with 100 iterations...
Best score at iteration 0: 725468.99; overall: 725468.99 (0s)
Best score at iteration 1: 804441.6; overall: 725468.99 (0s)
Best score at iteration 2: 824961.53; overall: 725468.99 (0s)
Best score at iteration 3: 898382.85; overall: 725468.99 (0s)
Best score at iteration 4: 897055.51; overall: 725468.99 (0s)
Best score at iteration 5: 763540.39; overall: 725468.99 (0s)
Best score at iteration 6: 754611.39; overall: 725468.99 (0s)
Best score at iteration 7: 740694.62; overall: 725468.99 (0s)
Best score at iteration 8: 739618.61; overall: 725468.99 (0s)
Best score at iteration 9: 746949.15; overall: 725468.99 (0s)
Best score at iteration 10: 740694.62; overall: 725468.99 (0s)
Best score at iteration 11: 750253.08; overall: 725468.99 (0s)
Best score at iteration 12: 727348.28; overall: 725468.99 (0s)
Best score at iteration 13: 731547.84; overall: 725468.99 (0s)
Best score at iteration 14: 727348.28; overall: 725468.99 (0s)
Best score at iteration 15: 718515.62; overall: 718515.62 (0s)
Best score at iteration 16: 727348.28; overall: 718515.62 (0s)
Best score at iteration 17: 719899.96; overall: 718515.62 (0s)
Best score at iteration 18: 684749.9; overall: 684749.9 (0s)
Best score at iteration 19: 719671.99; overall: 684749.9 (0s)
Best score at iteration 20: 684749.9; overall: 684749.9 (0s)
Best score at iteration 21: 684749.9; overall: 684749.9 (0s)
Best score at iteration 22: 684749.9; overall: 684749.9 (0s)
Best score at iteration 23: 684749.9; overall: 684749.9 (0s)
Best score at iteration 24: 668281.37; overall: 668281.37 (0s)
Best score at iteration 25: 668281.37; overall: 668281.37 (0s)
Best score at iteration 26: 668281.37; overall: 668281.37 (0s)
Best score at iteration 27: 668281.37; overall: 668281.37 (0s)
Best score at iteration 28: 668281.37; overall: 668281.37 (0s)
Best score at iteration 29: 668281.37; overall: 668281.37 (0s)
Best score at iteration 30: 668281.37; overall: 668281.37 (0s)
Best score at iteration 31: 668281.37; overall: 668281.37 (0s)
Best score at iteration 32: 668281.37; overall: 668281.37 (0s)
Best score at iteration 33: 668281.37; overall: 668281.37 (0s)
Best score at iteration 34: 668281.37; overall: 668281.37 (0s)
Best score at iteration 35: 668281.37; overall: 668281.37 (0s)
Best score at iteration 36: 668281.37; overall: 668281.37 (0s)
Best score at iteration 37: 668281.37; overall: 668281.37 (0s)
Best score at iteration 38: 668281.37; overall: 668281.37 (0s)
Best score at iteration 39: 668281.37; overall: 668281.37 (0s)
Best score at iteration 40: 668281.37; overall: 668281.37 (0s)
Best score at iteration 41: 668281.37; overall: 668281.37 (0s)
Best score at iteration 42: 668281.37; overall: 668281.37 (0s)
Best score at iteration 43: 668281.37; overall: 668281.37 (0s)
Best score at iteration 44: 668281.37; overall: 668281.37 (0s)
Best score at iteration 45: 668281.37; overall: 668281.37 (0s)
Best score at iteration 46: 668281.37; overall: 668281.37 (0s)
Best score at iteration 47: 668281.37; overall: 668281.37 (0s)
Best score at iteration 48: 668281.37; overall: 668281.37 (0s)
Stopping early due to 20 iterations of the same score.
ACO fitted.  Runtime: 0 minutes.  Best score: 668281.3692114921

Beginning ACO Optimization with 100 iterations...
Best score at iteration 0: 775390.1; overall: 775390.1 (0s)
Best score at iteration 1: 885962.49; overall: 775390.1 (0s)
Best score at iteration 2: 726420.18; overall: 726420.18 (0s)
Best score at iteration 3: 880384.95; overall: 726420.18 (0s)
Best score at iteration 4: 799468.6; overall: 726420.18 (0s)
Best score at iteration 5: 726420.18; overall: 726420.18 (0s)
Best score at iteration 6: 766951.4; overall: 726420.18 (0s)
Best score at iteration 7: 731997.72; overall: 726420.18 (0s)
Best score at iteration 8: 726283.37; overall: 726283.37 (0s)
Best score at iteration 9: 726420.18; overall: 726283.37 (0s)
Best score at iteration 10: 724611.44; overall: 724611.44 (0s)
Best score at iteration 11: 726420.18; overall: 724611.44 (0s)
Best score at iteration 12: 726283.37; overall: 724611.44 (0s)
Best score at iteration 13: 724611.44; overall: 724611.44 (0s)
Best score at iteration 14: 715136.57; overall: 715136.57 (0s)
Best score at iteration 15: 713699.13; overall: 713699.13 (0s)
Best score at iteration 16: 715136.57; overall: 713699.13 (0s)
Best score at iteration 17: 715136.57; overall: 713699.13 (0s)
Best score at iteration 18: 715136.57; overall: 713699.13 (0s)
Best score at iteration 19: 696726.59; overall: 696726.59 (0s)
Best score at iteration 20: 715136.57; overall: 696726.59 (0s)
Best score at iteration 21: 713587.74; overall: 696726.59 (0s)
Best score at iteration 22: 696726.59; overall: 696726.59 (0s)
Best score at iteration 23: 693705.55; overall: 693705.55 (0s)
Best score at iteration 24: 675295.57; overall: 675295.57 (0s)
Best score at iteration 25: 666529.07; overall: 666529.07 (0s)
Best score at iteration 26: 684939.04; overall: 666529.07 (0s)
Best score at iteration 27: 666529.07; overall: 666529.07 (0s)
Best score at iteration 28: 666529.07; overall: 666529.07 (0s)
Best score at iteration 29: 666529.07; overall: 666529.07 (0s)
Best score at iteration 30: 666529.07; overall: 666529.07 (0s)
Best score at iteration 31: 666529.07; overall: 666529.07 (0s)
Best score at iteration 32: 666529.07; overall: 666529.07 (0s)
Best score at iteration 33: 666392.26; overall: 666392.26 (0s)
Best score at iteration 34: 666392.26; overall: 666392.26 (0s)
Best score at iteration 35: 666529.07; overall: 666392.26 (0s)
Best score at iteration 36: 666392.26; overall: 666392.26 (0s)
Best score at iteration 37: 666392.26; overall: 666392.26 (0s)
Best score at iteration 38: 666392.26; overall: 666392.26 (0s)
Best score at iteration 39: 666392.26; overall: 666392.26 (0s)
Best score at iteration 40: 666392.26; overall: 666392.26 (0s)
Best score at iteration 41: 666392.26; overall: 666392.26 (0s)
Best score at iteration 42: 666392.26; overall: 666392.26 (0s)
Best score at iteration 43: 666392.26; overall: 666392.26 (0s)
Best score at iteration 44: 666392.26; overall: 666392.26 (0s)
Best score at iteration 45: 666392.26; overall: 666392.26 (0s)
Best score at iteration 46: 666392.26; overall: 666392.26 (0s)
Best score at iteration 47: 666392.26; overall: 666392.26 (0s)
Best score at iteration 48: 666392.26; overall: 666392.26 (0s)
Best score at iteration 49: 666392.26; overall: 666392.26 (0s)
Best score at iteration 50: 666392.26; overall: 666392.26 (0s)
Best score at iteration 51: 666392.26; overall: 666392.26 (0s)
Best score at iteration 52: 666392.26; overall: 666392.26 (0s)
Best score at iteration 53: 666392.26; overall: 666392.26 (0s)
Best score at iteration 54: 666392.26; overall: 666392.26 (0s)
Best score at iteration 55: 666392.26; overall: 666392.26 (0s)
Stopping early due to 20 iterations of the same score.
ACO fitted.  Runtime: 0 minutes.  Best score: 666392.259629674

Beginning ACO Optimization with 100 iterations...
Best score at iteration 0: 857230.94; overall: 857230.94 (0s)
Best score at iteration 1: 865477.55; overall: 857230.94 (0s)
Best score at iteration 2: 888071.88; overall: 857230.94 (0s)
Best score at iteration 3: 916103.29; overall: 857230.94 (0s)
Best score at iteration 4: 825349.44; overall: 825349.44 (0s)
Best score at iteration 5: 915775.09; overall: 825349.44 (0s)
Best score at iteration 6: 840176.01; overall: 825349.44 (0s)
Best score at iteration 7: 833182.93; overall: 825349.44 (0s)
Best score at iteration 8: 687450.69; overall: 687450.69 (0s)
Best score at iteration 9: 765471.01; overall: 687450.69 (0s)
Best score at iteration 10: 695595.78; overall: 687450.69 (0s)
Best score at iteration 11: 709547.32; overall: 687450.69 (0s)
Best score at iteration 12: 756948.6; overall: 687450.69 (0s)
Best score at iteration 13: 695342.3; overall: 687450.69 (0s)
Best score at iteration 14: 693615.5; overall: 687450.69 (0s)
Best score at iteration 15: 730712.44; overall: 687450.69 (0s)
Best score at iteration 16: 703433.44; overall: 687450.69 (0s)
Best score at iteration 17: 693478.7; overall: 687450.69 (0s)
Best score at iteration 18: 693615.5; overall: 687450.69 (0s)
Best score at iteration 19: 688284.18; overall: 687450.69 (0s)
Best score at iteration 20: 693615.5; overall: 687450.69 (0s)
Best score at iteration 21: 679161.08; overall: 679161.08 (0s)
Best score at iteration 22: 680124.51; overall: 679161.08 (0s)
Best score at iteration 23: 688804.55; overall: 679161.08 (0s)
Best score at iteration 24: 679161.08; overall: 679161.08 (0s)
Best score at iteration 25: 670394.58; overall: 670394.58 (0s)
Best score at iteration 26: 670394.58; overall: 670394.58 (0s)
Best score at iteration 27: 670257.77; overall: 670257.77 (0s)
Best score at iteration 28: 670257.77; overall: 670257.77 (0s)
Best score at iteration 29: 670257.77; overall: 670257.77 (0s)
Best score at iteration 30: 670257.77; overall: 670257.77 (0s)
Best score at iteration 31: 670257.77; overall: 670257.77 (0s)
Best score at iteration 32: 670257.77; overall: 670257.77 (0s)
Best score at iteration 33: 670257.77; overall: 670257.77 (0s)
Best score at iteration 34: 670257.77; overall: 670257.77 (0s)
Best score at iteration 35: 670257.77; overall: 670257.77 (0s)
Best score at iteration 36: 670257.77; overall: 670257.77 (0s)
Best score at iteration 37: 670257.77; overall: 670257.77 (0s)
Best score at iteration 38: 670257.77; overall: 670257.77 (0s)
Best score at iteration 39: 670257.77; overall: 670257.77 (0s)
Best score at iteration 40: 670257.77; overall: 670257.77 (0s)
Best score at iteration 41: 670257.77; overall: 670257.77 (0s)
Best score at iteration 42: 670257.77; overall: 670257.77 (0s)
Best score at iteration 43: 670257.77; overall: 670257.77 (0s)
Best score at iteration 44: 670257.77; overall: 670257.77 (0s)
Best score at iteration 45: 670257.77; overall: 670257.77 (0s)
Best score at iteration 46: 670257.77; overall: 670257.77 (0s)
Best score at iteration 47: 670257.77; overall: 670257.77 (0s)
Best score at iteration 48: 670257.77; overall: 670257.77 (0s)
Best score at iteration 49: 670257.77; overall: 670257.77 (0s)
Best score at iteration 50: 670257.77; overall: 670257.77 (0s)
Best score at iteration 51: 670257.77; overall: 670257.77 (0s)
Stopping early due to 20 iterations of the same score.
ACO fitted.  Runtime: 0 minutes.  Best score: 670257.7702044595

Calculamos el costo total del recorrido

def costo_total(ruta, costos):
    costo = 0
    for i in range(len(ruta) - 1):
        costo += MCTC[ruta[i]][ruta[i+1]]
    return costo

for i in range(len(MejorRuta)):
    print("El costo total de la ruta es: $", round(costo_total(MejorRuta[i], MCTC), 2), " COP para un total de ", Hormigas[i], "Hormigas")

El costo total de la ruta es: $ 729205.86  COP para un total de  10 Hormigas
El costo total de la ruta es: $ 668281.37  COP para un total de  20 Hormigas
El costo total de la ruta es: $ 668281.37  COP para un total de  80 Hormigas
El costo total de la ruta es: $ 666392.26  COP para un total de  100 Hormigas
El costo total de la ruta es: $ 670257.77  COP para un total de  120 Hormigas

#Cada vez que ejecutemos esta celda y el script en general borramos los archivos para que no se acumulen plot de otras seseiones
folder_path = "GifVCH"
for file_name in os.listdir(folder_path):
    file_path = os.path.join(folder_path, file_name)
    if os.path.isfile(file_path):
        os.remove(file_path)
# Lista para almacenar las rutas de las rutas de las imágenes PNG
ruta_imagenes = []

for i in range(len(MejorRuta)):
    texto = "La ruta es: "+"\n"  # texto inical para añadir las rutas al plot
    texto2 = "El costo total de la ruta es: $" + str(round(costo_total(MejorRuta[i], MCTC), 2))+"COP"  # Texto para añadir el valor de la ruta

    current_sol = MejorRuta[i]  # Obtenbemos la soulcion por iteracion
    for j in range(len(current_sol)):
        texto += str(ciudades[current_sol[j]])+"\n"

    # Aqui obtenemos la ruta segun sea la solucion
    rutacol = [df.iloc[i]['geometry']for i in MejorRuta[i]]
    #creamos el plot de la ruta
    fig, ax = plt.subplots(figsize=(14, 14))
    mapa.plot(ax=ax, color='lime', edgecolor='black')
    df.plot(ax=ax, markersize=50, color='red')
    for i, txt in enumerate(df['Ciudad']):
        ax.annotate(txt, (df['geometry'][i].x,
                    df['geometry'][i].y), fontsize=12)
    df.plot(ax=plt.gca(), color='red', markersize=50, zorder=3)
    plt.plot([p.x for p in rutacol], [p.y for p in rutacol],
             color='blue', linestyle='--', linewidth=2, dashes=[2, 2], zorder=2)
    plt.text(250000, 1000000, texto, ha='center', va='center', size=14,
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.9))
    plt.text(550000, 2000000, texto2, ha='center', va='center', size=14,
             bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.9))
    plt.title("Problema del viajero: Colonias de Hormigas")
    # Generar un nombre de archivo único con un número aleatorio
    nombre_archivo = f"GifVCH/imagen_{i}_{random.randint(0, 100000)}.png"
    texto = ""
    # Agregar la ruta de la imagen PNG a la lista
    ruta_imagenes.append(nombre_archivo)

    # Guardar la figura como imagen PNG
    plt.savefig(nombre_archivo)

# Crear archivo GIF a partir de las imágenes PNG
with imageio.get_writer('GifVCH/animacionCH.gif', mode='I', duration=0.8) as writer:
    for ruta in ruta_imagenes:
        image = imageio.imread(ruta)
        writer.append_data(image)

Problema del Viajero Colombian Version: Algoritmos geneticos

#Variables necesarias para la creacion de gif
num_imagenes = 10
imagenes = []

def tsp_fitness(solution, s_idx):
    #Realizamos unos ajustes a soltion puesto que no se encuentran en formato que nos deje iterar
    solution = solution.tolist()
    solution = [int(x) for x in solution]

    # Verificar que no haya ciudades repetidas en la solución
    if len(set(solution)) != len(solution):
        return 0  # retorna aptitud cero si hay ciudades repetidas en la solución

    # Verificar que la solución contenga todas las ciudades
    if set(solution) != set(range(len(ciudades))):
        return 0  # retorna aptitud cero si falta alguna ciudad en la solución

    # Calcular el costo total del recorrido
    distance = 0
    for i in range(len(solution)-1):
        distance += MCTC[solution[i]][solution[i+1]]
    distance += MCTC[solution[-1]][solution[0]]  # vuevle al punto de inicio

    # Retornar la inversa de la distancia como aptitud (para maximizar)
    return 1/distance

Generaciones = [100, 150, 270, 500, 1120, 1500]
mejor_soluciones = []
sol_per_pop = 15
soluciones_iniciales = [list(set(np.random.permutation(idx_ciudades)))
                        for _ in range(sol_per_pop)]
for i in range(len(Generaciones)):
    num_generations = Generaciones[i]
    # Definir el modelo de optimización genética con PyGAD
    ViajeroGA = pygad.GA(num_generations=num_generations,
                        num_parents_mating=4,
                        sol_per_pop=sol_per_pop,
                        num_genes=len(ciudades),
                        initial_population=soluciones_iniciales,
                        mutation_percent_genes=30,
                        fitness_func=tsp_fitness,
                        mutation_type="scramble",
                        init_range_low=1,
                        init_range_high=15,
                        parent_selection_type="rank")
    # Ejecutar la optimización genética
    ViajeroGA.run()

    # Obtener la mejor solución encontrada
    solution, solution_fitness, solution_idx = ViajeroGA.best_solution()
    solution_indices = [int(x) for x in solution.tolist()]
    mejor_soluciones.append(solution_indices)
    MejorCiudades = []
    for i in solution_indices:
        MejorCiudades.append(ciudades[i])

    # Imprimir la mejor solución encontrada
    print("Mejor solución encontrada en:",
          num_generations, "iteraciones es:", MejorCiudades)

    

Mejor solución encontrada en: 100 iteraciones es: ['Cúcuta', 'Bucaramanga', 'Bogota', 'Armenia', 'Tuluá', 'Pasto', 'Palmira', 'Pereira', 'Manizales', 'Medellín', 'Montería', 'Cartagena', 'Soledad', 'Barranquilla', 'Valledupar']
Mejor solución encontrada en: 150 iteraciones es: ['Montería', 'Bogota', 'Manizales', 'Medellín', 'Tuluá', 'Palmira', 'Pasto', 'Pereira', 'Armenia', 'Valledupar', 'Barranquilla', 'Soledad', 'Cartagena', 'Bucaramanga', 'Cúcuta']
Mejor solución encontrada en: 270 iteraciones es: ['Manizales', 'Armenia', 'Pereira', 'Tuluá', 'Palmira', 'Pasto', 'Bogota', 'Bucaramanga', 'Cúcuta', 'Valledupar', 'Cartagena', 'Barranquilla', 'Soledad', 'Montería', 'Medellín']
Mejor solución encontrada en: 500 iteraciones es: ['Bogota', 'Armenia', 'Tuluá', 'Palmira', 'Pasto', 'Pereira', 'Manizales', 'Medellín', 'Montería', 'Cúcuta', 'Cartagena', 'Soledad', 'Barranquilla', 'Valledupar', 'Bucaramanga']
Mejor solución encontrada en: 1120 iteraciones es: ['Bogota', 'Armenia', 'Palmira', 'Pasto', 'Tuluá', 'Pereira', 'Medellín', 'Montería', 'Soledad', 'Barranquilla', 'Cartagena', 'Valledupar', 'Cúcuta', 'Bucaramanga', 'Manizales']
Mejor solución encontrada en: 1500 iteraciones es: ['Valledupar', 'Bogota', 'Armenia', 'Pasto', 'Palmira', 'Tuluá', 'Pereira', 'Manizales', 'Medellín', 'Montería', 'Barranquilla', 'Soledad', 'Cartagena', 'Cúcuta', 'Bucaramanga']

Calculamos el costo total del recorrido

def costo_total(ruta, costos):
    costo = 0
    for i in range(len(ruta) - 1):
        costo += MCTC[ruta[i]][ruta[i+1]]
    return costo

for i in range(len(mejor_soluciones)):
    print("El costo total de la ruta es: $",round(costo_total(mejor_soluciones[i],MCTC),2)," COP para un total de ", Generaciones[i], "Iteraciones")

El costo total de la ruta es: $ 648272.0  COP para un total de  100 Iteraciones
El costo total de la ruta es: $ 749751.41  COP para un total de  150 Iteraciones
El costo total de la ruta es: $ 653715.53  COP para un total de  270 Iteraciones
El costo total de la ruta es: $ 655548.19  COP para un total de  500 Iteraciones
El costo total de la ruta es: $ 641403.75  COP para un total de  1120 Iteraciones
El costo total de la ruta es: $ 662283.01  COP para un total de  1500 Iteraciones

#Cada vez que ejecutemos esta celda y el script en general borramos los archivos para que no se acumulen plot de otras seseiones
folder_path = "GifVGA"
for file_name in os.listdir(folder_path):
    file_path = os.path.join(folder_path, file_name)
    if os.path.isfile(file_path):
        os.remove(file_path)
# Lista para almacenar las rutas de las rutas de las imágenes PNG
ruta_imagenes = []

for i in range(len(mejor_soluciones)):
    texto = "La ruta es: "+"\n" # texto inical para añadir las rutas al plot
    texto2 = "El costo total de la ruta es: $"+str(round(costo_total(mejor_soluciones[i], MCTC), 2))+"COP" #Texto para añadir el valor de la ruta

    current_sol = mejor_soluciones[i] #Obtenbemos la soulcion por iteracion
    for j in range(len(current_sol)):
        texto += str(ciudades[current_sol[j]])+"\n"

    rutacol = [df.iloc[i]['geometry']for i in mejor_soluciones[i]] #Aqui obtenemos la ruta segun sea la solucion
    #creamos el plot de la ruta
    fig, ax = plt.subplots(figsize=(14, 14))
    mapa.plot(ax=ax, color='lime', edgecolor='black')
    df.plot(ax=ax, markersize=50, color='red')
    for i, txt in enumerate(df['Ciudad']):
        ax.annotate(txt, (df['geometry'][i].x,
                    df['geometry'][i].y), fontsize=12)
    df.plot(ax=plt.gca(), color='red', markersize=50, zorder=3)
    plt.plot([p.x for p in rutacol], [p.y for p in rutacol],
             color='blue', linestyle='--', linewidth=2, dashes=[2, 2], zorder=2)
    plt.text(250000, 1000000, texto, ha='center', va='center', size=14,
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.9))
    plt.text(550000, 2000000, texto2, ha='center', va='center', size=14,
             bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.9))
    plt.title("Problema del viajero: Colonias de Hormigas")

    # Generar un nombre de archivo único con un número aleatorio
    nombre_archivo = f"GifVGA/imagen_{i}_{random.randint(0, 100000)}.png"
    texto=""
    # Agregar la ruta de la imagen PNG a la lista
    ruta_imagenes.append(nombre_archivo)

    # Guardar la figura como imagen PNG
    plt.savefig(nombre_archivo)

# Crear archivo GIF a partir de las imágenes PNG
with imageio.get_writer('GifVGA/animacionGA.gif', mode='I', duration=0.8) as writer:
    for ruta in ruta_imagenes:
        image = imageio.imread(ruta)
        writer.append_data(image)

Animaciones de las rutas

from IPython.display import Image, display

gif1 = Image(filename='GifVCH/animacionCH.gif')
gif2 = Image(filename='GifVGA/animacionGA.gif')
display(gif1, gif2)

<IPython.core.display.Image object>

<IPython.core.display.Image object>

Conclusiones del problema del viajero

La optimizacion del problema por Colonias de hormigas nos da las rutas mas baratas, tambien se nota que por este metodo convergen mas rapido las soluciones a rutas Optimas

Referencias

Repositorio del proyecto: https://github.com/oyeu/Trabajo1-Optimizaci-n-heuristica